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SECTION  1 
INTRODUCTION 


The  WEDCOM  code  is  a digital  computer  program  that  provides 
calculations  of  ELF,  VLF,  and  LF  electric  and  magnetic  field  strengths 
in  ambient  and  nuclear  disturbed  environments.  The  code  is  intended 
for  use  when  a relatively  detailed  analysis  of  the  propagation  between 
two  terminals  is  required  in  nuclear  disturbed  environments.  The  code 
can  be  used  alone  to  study  effects  of  weapon,  atmosphere,  and  propaga- 
tion parameters  on  received  Signals  or  can  be  used  in  conjunction  with 
receiver  antenna  and  signal  processing  models  to  evaluate  system 
performance. 

Operational  information  for  WEDCOM  IV  is  given  in  Volume  1 
(User's  Manual).  In  this  volume  a description  of  computer  routines  is 
given  for  use  in  conjunction  with  FORTRAN  source  listings.  A mathe- 
matical description  of  computational  models  used  in  the  routines  is 
given  in  Volume  3. 

The  WEDCOM  IV  code  consists  of  a main  driver  routine  and  a 
number  of  subroutines  and  functions.  The  FORTRAN  names  and  a brief 
description  of  the  purpose  of  each  routine  are  given  in  Table  1.  As 
indicated  in  Table  1,  many  of  the  environment  and  geometry  routines 
are  taken  directly  from  the  WEPH  VI  code.  These  routines  are  des- 
cribed in  Reference  1 and  are  not  further  discussed  here. 

Section  2 of  this  volume  contains  intermediate  level  descrip- 
tions of  each  routine.  Except  for  the  main  driver  routine,  which  is 
described  first,  the  descriptions  are  presented  in  alphabetical  order. 
A number  of  the  routines  are  the  same  as  routines  used  in  a code  de- 
veloped for  the  MEECN  office  designated  DCOM  (Reference  2)  and  the 


descriptions  given  in  Section  2 for  these  routines  have  been  taken 
from  the  DCOM  documentation.  Table  2 lists  WEDCOM  IV  routines, 
library  routines,  and  labeled  common  blocks  used  by  each  routine. 


Table  1.  Description  of  WEDCOM  IV  routines. 


Name  of  Routine 


Control 


ABFACR 


>ANLYT2 


APPROX 


ARRLIM 


ATMOSF 


ATMOSU 


AZELR 


*BFIELD 


BHRHS 


Description 

Determines  sequence  of  calculations 

Computes  terms  used  in  an  approximate  mode 
equation 

Computes  the  Airy  functions  and  their  deriva- 
tives 

Computes  solution  to  differential  equation 
used  in  deionization  above  100  km 

Computes  eigenvalue  solutions  to  the  approximate 
mode  equation  by  employing  the  Newton  iteration 
algorithm  after  an  appropriate  initial  value  is 
derived 

Determines  debris  location  and  size  from 
debris  marker  particle  locations 

Provides  properties  of  the  normal  or  heaved 
atmosphere  above  100  km 

Provides  properties  of  the  normal  atmosphere 

Computes  true  azimuth,  elevation,  and  slant 
range  of  one  point  relative  to  another  point 
(points  specified  by  vectors) 

Determines  azimuth  of  vector  from  the  x-axis 
measured  in  x-y  plane 

Computes  properties  of  magnetic  dipole  field 
at  a point 

Computes  starting  and  stopping  altitudes  for 
the  reflection  coefficient  integration  for  LF 
propagation 


* Unmodified  WEPH  VI  routine. 


(Continued) 
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Table  1.  (Continued) 


Name  of  Routine 

* BLKATM 

* BLKCHM 

* BLKDEP 

* BLKUV 

BUDDEN 

CANG 

*C11EMI)2 

Cl  II. MIX) 

♦CHEMEF 

*CIIEMQ 

♦CHMION 

♦CHXLOS 

♦CHXSPC 

♦CLOSE 

♦CMFEDT 

CON FAC 
CONMAP 


Description 

Block  data  for  atmospheric  routines 

Block  data  for  E-  and  F-region  chemistry 
rout i nes 

Block  data  for  heavy  particle  energy  deposi- 
tion calculations 

Block  data  for  UV  energy  deposition 
calculations 

Driver  routine  for  determining  value  of  reflec- 
tion coefficient  for  ELF  or  VLF  propagation 

Calculates  the  argument  of  a complex  number 

Driver  routine  for  D-region  chemistry  module 


Driver  routine  for  steady-state  D-region 
chemistry 

Determines  chemistry  and  ionization  for  E-  and 
F-regions 

Determines  steady-state  ionization  above  100  km 

Determines  species  concentrations  in  E-  and  F- 
regions 

Determines  energy  emitted  in  charge  exchange 
particles 

Determines  energy  histogram  for  charge  exchange 
particles 

Determines  location  of  closest  point  of  approach 
between  two  vectors 

Edits  burst  list  for  bursts  to  be  used  in  E-  and 
F-region  chemistry  calculations 

Computes  ionospheric  convergence  factor 

Contains  a modified  version  of  the  ITS  world 
conductivity  map 


Unmodified  WEPH  VI  routine. 


(Continued) 


Table  1.  (Continued) 


Name  of  Routine 


Description 


*CONSET 


Sets  constants  used  in  low-altitude  fireball/ 
debris  models 


*CONSPC 

‘COORDV 

*COORDX 


‘CROSS 

‘DATE 

‘DBMHT 

‘DBMLT 


Determines  energy  histogram  for  loss-cone 
particles 

Determines  vector  coordinates  from  latitude, 
longitude,  and  altitude 

Transforms  coordinates  of  a point  in  one 
coordinate  system  to  any  of  the  other  three 
systems;  the  coordinate  systems  are 

1.  Ground  range,  azimuth,  altitude 

2.  X (east),  Y (north),  Z (vertical) 

3.  Slant  range,  azimuth,  elevation 

4.  Elevation,  azimuth,  altitude 


Calculates  cross  product  for  two  vectors 

Computes  Gregorian  date  and  zone  time  after 
specified  time  interval 

Determines  time-dependent  debris  region  quanti- 
ties for  high-altitude  bursts 


Determines  time-dependent  debris  region  quanti- 
ties for  low-altitude  bursts  (RANC  model) 


‘DBMLTR 


Determines  time-dependent  debris  region  quanti- 
ties for  low-altitude  bursts  (ROSCOE  model) 


‘DBSTAM 

‘DEBRIS 

‘DF.DEP 

DELTAT 

'DISTF 

‘DOT 


Driver  routine  for  post-stabilization  debris 
geometry 

Driver  routine  to  determine  energy  deposition 
by  heavy  particles 

Determines  delayed  energy  deposition  rate 

Computes  correction  to  an  approximate  root  of 
the  mode  equation 

Computes  normalized  debris  radial  distribution 
parameter 

Calculates  vector  dot  product 


* Unmodified  WEPH  VI  routine. 


(Continued) 
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Table  1.  (Continued) 


Name  of  Routine 


*DRATE 


‘DTNEP 


*DTNEQ 


IHISTMI 


‘EDEPB 


‘EDEPC 


*EDEPG 


‘EDF.PM 


►EDEPNB 


‘EDF.PND 


FIGI.NV 


'EIGSRH 


Description 


Determines  reaction  rate  coefficients  for  D- 
region  chemistry  solutions 

Computes  positive  ion  and  electron  concentra- 
tions at  a point  at  or  below  100  km 

Computes  D-region  positive  ion  and  electron 
concentrations  due  to  prompt  radiation 

Computes  steady-state  D-region  positive  ion 
and  electron  concentrations  due  to  delayed 
radiation 

Computes  dust  scaling  quantities 

Determines  energy  deposition  rate  due  to  beta 
particles 

Determines  energy  deposition  rate  due  to 
Compton  electrons 

Determines  energy  deposition  rate  due  to  gamma 
rays 

Determines  energy  deposition  rate  due  to 
Bremsst rah  lung  radiation 

Determines  energy  deposition  rate  due  to  neutron 
decay  beta  particles 

Determines  energy  deposition  rate  due  to  neutron 
elastic  collisions  and  capture  reactions 

Computes  vertical  and  horizontal  field  strength 
for  Jth  hop  in  IF  propagation 

Calculates  the  exponential  integral  function 
for  positive  arguments 

Calculates  the  eigenvalue  and  derivative  to  the 
mode  equation  using  Newton-Raphson  procedure 

Computes  VLF  mode  characteristics  from  NELC 
mode  search  model 


* Unmodified  WEPH  VI  routine. 


(Continued) 
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Table  1.  (Continued) 


Name  of  Routine 


Description 


F.LCOL 

elde.mj 


r.LFNion 


Computes  colatitude  and  east  longitude  of  a 
point  on  the  earth 

Driver  routine  for  chemistry  and  ionization 
cal  cut  at  ions 

Determines  elevation  of  vector  from  the  x-y 
plane 

Driver  routine  for  the  calculation  of  ELF 


EMATRX  Computes  elements  of  a susceptibility  matrix 

including  magnetic  field  effects  and  earth 
curvature  for  ELF  propagation. 

*ENEF  Driver  routine  for  E-  and  F-region  chemistry 

module 


ENV1RM 

ENNU 

EPRIME 
*EQRVT 
ER I NT 

ESMAT 

F.START 

ETGAIN 


Driver  routine  for  calculating  environment 
quantities  along  a given  path 

Computes  electron  and  ion  densities  and  colli- 
sion frequencies  through  an  exponential  inter- 
polation between  data  points 

Computes  reflection  coefficient  derivatives 
with  respect  to  altitude 

Determines  effective  ambient  photodissociation 
rates  for  D- region 

Finds  reflection  coefficients  using  Runge-Kutta 
integration  for  ELF  propagation 

Computes  coefficient  matrix  for  the  reflection 
coefficient  differential  equation 

Computes  starting  values  for  reflection  coeffi- 
cient integration  for  ELF  propagation 

Calculates  height  gain  function  for  ELF 


FXCFAC 

EUXF1T 


Calculates  anisotropic  excitation  factors 

Determines  total  UV  radiation  and  fraction 
radiated  in  five  radiation  groups 


* Unmodified  WEPH  VI  routine. 


(Continued) 
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Table  1.  (Continued) 


Name  of  Routine 

Description 

*EX1NTP 

Provides  exponential  interpolation 

*E2 

Evaluates  function  used  in  beta  particle 
energy  deposition  routine 

*EBMHI 

Determines  t ime- independent  fireball  quanti- 
ties for  high-altitude  bursts 

*FBN1HT 

Determines  time-dependent  fireball  quanti- 
ties for  high-altitude  bursts 

*FBML1 

Determines  t ime- independent  fireball  quanti- 
ties for  low-altitude  bursts  (RANG  model) 

*FBMLT 

Determines  time-dependent  fireball  quanti- 
ties for  low-altitude  bursts  (RANG  model) 

*FHMNMX 

Determines  top  and  bottom  altitude  of  a tube 
fireball  and  altitude  of  fireball  at  magnetic 
equator  if  fireball  reaches  equator 

12 
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Table  1.  (Continued) 

Name  of  Routine 

GUESS 

*HDPART 
*HEAVE 

*HELP 
*HPCHEM 
HTGAIN 

*HTOS 
*HYDMRG 

* I N1TAL 
INPUT 

* 10NLEK 

IONOSU 
ISOREF 

1S0RLF 

*JULIAN 

*LAGRAN 

* Unmodified  WEPH  VI  routine. 


Description 

Computes  eigenvalue  assuming  an  isotropic 
ionosphere 

Driver  routine  to  partition  hydro  energy 
into  heavy  particle  energy  sources 

Computes  .the  vertical  velocity  and  time  his- 
tory of  the  altitude,  expansion  ratio,  and 
density  scale  height  of  a cell  in  a vertically 
heaving  atmosphere 

Prints  an  error  message  and  can  continue  or 
terminate  the  problem 

Partitions  energy  lost  by  heavy  particles  and 
determines  change  in  neutral  and  ion  species 

Calculates  antenna  height  gain  factors  for 
VLF  propagation 

Calculates  slant  ranges  along  vector  to  points 
that  are  at  specified  altitude 

Determines  t ime- independent  parameters  for 
fireball  formed  by  hydro  merge 

Determines  D-region  neutral  species  concen- 
trations at  end  of  Phase  2 following  a burst 

Reads  and  writes  input  quantities  and  performs 
preliminary  geometry  calculations 

Determines  energy  in  ion  leak  particles  and 
spatial  distribution  parameter 

Provides  normal  ionospheric  quantities 

Computes  VLF  isotropic  ionospheric  reflection 
coefficient 

Computes  LF  isotropic  ionospheric  reflection 
coefficient 

Computes  Gregorian  calendar  date  to  a Julian 
day  number 

Computes  the  initial  location  and  final  prop- 
erties of  a Lagrangian  cell 


(Continued) 
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Table  1.  (Continued) 


Name  of  Routine 
*LEKSPC 

LENTRP 

LFHOP 

LINKER 

*LOCLAX 
* LOS CON 

*MAGFIT 

*MAGFLD 

MATINV 

*MATMUL 

MDHNKL 

NftATRX 

MODCON 

‘MODEL 

‘MODELT 

* Unmodified  WEPH 


Description 

Determines  energy  histogram  for  ion  leak 
particles 

Interpolates  linearly  to  obtain  electron  and 
positive  ion  density  values  for  normal  and 
disturbed  conditions 

Driver  routine  for  the  calculation  of  LF  propa- 
gation 

Used  on  the  Honeywell  machine  for  linking 
(overlay) 

Calculates  vector  transformation  matrix 

Determines  energy  in  loss  cone  particles  and 
spatial  distribution  parameter 

Fits  a dipole  magnetic  field  to  the  local 
magnetic  field  at  specified  point 

Determines  location  of  the  point  on  a magnetic 
field  line  that  is  at  a specified  altitude 

Solves  a set  of  linear  simultaneous  equations 

Performs  matrix  multiplication  between  two 
real  matrices 

Computes  terms  used  in  forming  modified 
Hankel  functions 


Computes  elements  of  a susceptibility  matrix 
including  magnetic  field  effects  and  earth 
curvature  for  VLF  and  LF  propagation 

Computes  mode  conversion  coefficients 

Driver  routine  for  ROSCOE  low-altitude  fire- 
bal 1/debris  model 

Determines  time-dependent  fireball/debris 
quantities  for  low-altitude  bursts  (ROSCOE 
model ) 


VI  routine. 


(Continued) 
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Table  1.  (Continued) 

Name  of  Routine 


MODEZ 

*MODLON 


‘OFFSET 

‘ONF.MGS 
OUT  I ON 


*PCHEM 

*PEDEP 

*PHEAVE 


Description 

Computes  VLF  mode  characteristics  from  TEMPO 
mode  search  model 

Determines  t ime- independent  f i reba 1 1 /debr i s 
quantities  for  low-altitude  bursts  (ROSCOE 
model ) 

Determines  effects  of  shock  wave  on  low- 
altitude  firebal Is 

Calculates  magnetic  field  at  specified  point 

Computes  the  ionospheric  profile  index  of 
refraction  values,  the  reference  altitude  for 
VLF  reflection,  and  the  top  and  bottom  alti- 
tudes for  VLF  and  VLF  reflection  coefficient 
calculations 

Determines  change  in  species  densities  at 
altitudes  above  100  km  due  to  prompt  radiation 

Computes  the  energy  deposited  due  to  prompt 
radiation 

Determines  the  heave  disturbed  atmosphere 
prof i le 


PHENOM 


Driver  routine  for  the  RANC  phenomenology 
(fireball,  debris,  and  dust)  model  module 


*PHEN0R  Driver  routine  for  the  ROSCOE  phenomenology 

(fireball,  debris,  and  dust)  model  module 

‘PHOTOD  Computes  the  fraction  of  Oj  remaining  after 

thermal  radiation 


*PHOTOR  Calculates  the  negative  ion  photodetachment 

and  photodissociation  rates  due  to  thermal 
radiation 


‘PIONF 

‘PLINTP 


Obtains  the  initial  ion  concentrations  at  a 
point  above  100  km 

Obtains  the  power  law  interpolation, 
y = axb 


* Unmodified  WEPH  VI  routine. 


(Continued) 
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Table  1.  (Continued) 

Name  of  Routine  Description 


*PMASS 

*PSDM 

QSUBI 

RF.FLCO 

*RADMRG 

*RADOUT 

*RATE 

RF.ORDR 

RGND 

* R 1 CATT 

RINT 

*RI NTER 

ROOT 

RPR  IMF! 

RSTART 


Computes  the  mass  penetrated  between  two 
points  assuming  that  the  air  density  varies 
exponentially  with  altitude 

Driver  routine  for  post-stabilization  debris 
module  for  specified  debris  region 

Computes  the  normalized  ionospheric  surface 
impedance 

Driver  routine  for  determining  value  of  re- 
flection coefficient  for  LF  propagation 

Determines  t ime- independent  quantities  for 
fireball  formed  by  radiation  merge 

Computes  the  radiation  power  for  a burst  of 
a given  yield  at  a specified  altitude 

Calculates  the  reaction  rate  coefficient  for 
a specified  reaction 

Reorders  environmental  data  for  use  by  propa- 
gation modules 

Calculates  the  ground  reflection  coefficients 
for  ELF  and  VLF  propagation 

Computes  the  transient  electron  density 
based  on  an  approximate  solution  to  the 
Ricatti  equation 

Calculates  reflection  coefficients  using 
Runge-Kutta  integration  for  VLF  and  LF 
propagation 

Computes  intersections  between  ray  path  and 
right  circular  and  skewed  spheroids,  cones, 
magnetically  defined  tubes,  and  torus 

Computes  root  of  approximate  characteristic 
equation,  using  Newton's  procedure 

Computes  reflection  coefficient  derivative 
with  respect  to  altitude 

Computes  starting  values  for  reflection  coeffi- 
cient integration  for  VLF  and  LF  propagation 


* Unmodified  WEPH  VI  routine. 
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Table  1.  (Continued) 

Name  of  Routine 

RUNG 

*SEPA 

SMAT 

‘SOIL 

‘SOLCYC 

‘SOLORB 

‘SOLVE 

*SOLZEN 

‘SPCM1N 

‘SPECDP 

‘SPECDQ 

*STRIAT 

‘SUBVEC 

*SUB2 

*SUB3 

*SUB4 

*SUB5 

*SUB6 


Description 

Performs  Runge-Kutta  integration  in  ground- 
wave  calculation 

Calculates  the  angle  between  two  vectors 

Computes  coefficient  matrix  for  the  reflection 
coefficient  differential  equation 

Provides  the  mass  density  and  dielectric 
constant  for  a given  soil 

Computes  the  10.7-cm  solar  flux 

Computes  the  latitude  and  longitude  of  the 
subsolar  point 

Solves  a set  of  linear  algebraic  equations 
using  the  Gauss-Jordan  method 

Calculates  the  cosine  of  the  solar  zenith 
angle 

Computes  the  minor  neutral  atmospheric  species 

Computes  the  neutral  species  concentrations 
at  the  end  of  Phase  3 

Computes  the  modification  of  the  neutral 
species  due  to  delayed  radiation 

Provides  the  striation  ionization  parameters 
in  a late-time,  high-altitude  fireball 

Calculates  the  difference  between  two  vectors 

Determines  region  center  altitude  for  ROSCOE 
fireball  model 

Determines  debris  parameters  for  ROSCOE  model 

Determines  smoothing  parameter  for  ROSCOE  low- 
altitude  fireball  model 

Determines  region  altitude  for  ROSCOE  fireball 
model 

Determines  region  radii  for  ROSCOE  fireball 
model 


* Unmodified  WEPH  VI  routine. 
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Table  1.  (Continued) 

Name  of  Routine 


Description 


*SIJB10 

*SUB1 1 

*SUB1 2 

♦SIJB1 3 
*SUB14 
SURFQ 

TEM 

*TEXK 

TF1F.LD 

♦TOROID 

♦TUBE 

UANTD 

♦UNITV 

♦UVWAV 

♦UVWDH 

♦UVWDL1 

* Unmodified 


Determines  time  for  region  radius 
specified  factor  during  power-law 


to  grow  by 
expansion 


Determines  time  for  region  radius  to  grow  by 
specified  factor  during  deceleration  phase 


Determines  time  when  radius  reaches  specified 
value  in  a specified  direction 


Determines  vortex  radii  at  specified  time 


Determines  offset  of  fireball 


Computes  normalized  surface  impedance  for 
vertical  or  horizontal  polarization 

Calculates  transverse  electromagnetic  mode 
field  value  at  reveiver  for  ELF  propagation 

Computes  the  excitation  temperature  for  a 
given  energy  density  equilibrated  between  the 
N.  (VIB),  Of'D),  N(?l>),  and  the  electron  kinetic 
energy 

Computes  total  electric  field  from  hop  fields 
for  LF  propagation 

Computes  the  intersection  points  of  a ray 
path  with  a toroidal  region 

Computes  the  intersection  point  of  a ray  path 
with  a magnetically  contained  region 

User  supplied  routine  to  determine  transmitter 
antenna  orientation  as  a function  of  time 


Computes  a unit  vector  from  a given  input 
vector 

Determines  average  wind  velocity  over  a verti- 
cal mixing  length 

Determines  horizontal  wind  field  above  25-km 
altitude 

Determines  average  horizontal  wind  field  be- 
low 25-km  altitude  for  months  of  January, 

February,  March,  and  April 

WEPH  VI  routine.  (Continued) 
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Table  1.  (.Continued) 

Name  of  Routine 


— 


*UVWLD2 

*UVWLD3 

*U25 

*VECLIN 

*VECM 
* VECSUM 


Description 

Determines  average  horizontal  wind  field  be- 
low 25-km  altitude  for  months  of  May,  June, 
July,  and  August. 

Determines  average  horizontal  wind  field  be- 
low 25-km  altitude  for  months  of  September, 
October,  November,  and  December 

Determines  average  horizontal  wind  field  at 
25-km  altitude 

Computes  the  linear  combination  of  two  vectors 
each  multiplied  by  an  arbitrary  constant 

Computes  the  product  of  a vector  with  a scalar 

Computes  the  sum  of  two  arbitrary  vectors  and 
adds  it  to  a given  vector 


VFUNC1  Computes  a function  used  in  ROSCOE  fireball 

model 

*VFUNC2  Computes  a function  used  in  ROSCOE  fireball 

model 


VLFMCV  Determines  VLF  electric  and  magnetic  field 

strengths  including  effects  of  mode  conversion 


VLFMOD 

Driver  routine  for  calculation  of  VLF 
propagation 

VLI'WKB 

Determines  VLF  electric  and  magnetic  field 
strengths  neglecting  mode  conversion 

VOLUME 

Computes  the  volume  of  a given  geometrical 
region 

WDNATN 

Driver  routine  that  determines  ambient  atmos- 
phere and  ionosphere  properties 

WFCTVL 

Determine  F function  values  at  mesh  points 
for  NOSC  mode  search  model 

WFDFDT 

Determines  F function  values  for  NOSC  code 
search  model 

WF1NDF 

Used  in  finding  roots  in  NOSC  mode  search  mod' 

* Unmodified  WF.PH  VI  routine.  (Continued) 
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Table  1.  (Continued) 

Name  of  Routine 

WFSINT 

WFZERO 

WIN1FS 

WTNILG 

WINIRB 

WLGDER 

WLGRNG 

WNOMES 

*1V0BI) 

*WOG!) 

*1V0GP 
*WOG  1 
*U'ONI) 

*KONP 
*WON  1 
*NOXC 

* Unmodified  WLI’II  VI 


Description 

Performs  integration  for  the  ionisphere  re- 
flection matrix  for  NOSC  mode  search  model 

Determines  zero's  of  complex  function  for 
NOSC  mode  search  model 

Multiple  entry  of  routine  WFSINT 

Multiple  entry  of  routine  WLGRNG 

Multiple  entry  of  routine  WRBARS 

Multiple  entry  of  routine  WLGRNG 

Determines  Lagrange  interpolation  of  reflec- 
tion coefficients  for  NOSC  mode  search  model 

Determines  location  of  zeros  of  the  function 
F for  NOSC  mode  search  routine 

Computes  the  beta  energy  deposition  coeffi- 
cient, particle  release  rate,  and  average 
energy 

Computes  the  delayed  gamma-rav  energy 
deposition  parameters 

Provides  prompt  gamma  energy  deposition 
parameters 

Initializes  the  gamma -ray  energy  deposition 
parameters 

Provides  the  time-dependent  neutron  energy 
deposition  parameters 

Provides  prompt  neutron  energy  deposition 
parameters 

Initializes  the  neutron  energy  deposition 
parameters 

Computes  the  X-ray  energy  containment 


rout i ne. 


(Continued) 
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Table  1.  (Continued) 

Name  of  Routine  Description 

*WOXP  Computes  the  prompt  X-ray  energy  deposition 

parameters 

*W0X1  Initialises  the  neutron  energy  deposition 

parameters 

WQIIAD  Determines  roots  of  a quadratic  equation  for 

NOSC  mode  search  routine 

WRBARS  Determines  values  for  variables  used  in  RBAR 

matric  for  NOSC  mode  search  program 

KSfiT-Ci  Multiple  entry  of  routine  WLGRNG 

KSCTkll  Determines  given  points  in  Lagrange  interpo- 

lation for  NOSC  mode  search  program 

Xl-T.R  Transfers  one  array  to  another 

*XMAG  Computes  the  absolute  magnitude  (length)  of 

a vector 

*ITTOUT  Converts  a local  time  and  Gregorian  calendar 

datc  to  the  local  time  and  Gregorian  calen- 
dar date  to  Greenwich 


t, 


Table  2.  WEDCOM  IV  routines,  library  routines,  and  common  blocks 
called  directly  in  WEDCOM  IV  routines. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

Control 

AIRY 

FXOPT* 

ATMOSN 

ATMOSU 

LENTRY** 

ATMOUP 

AZELR 

SQRT 

BREG 

COORDV 

BURST 

COORDX 

CINPUT 

LINKER 

CMPLXJ 

MA6FIT 

CONBB 

VOLUME 

CONST 

ZTTOUT 

DBREG 

DCREG 

DEPDAT 

DEVICE 

EDTOVL 

FBREG 

FDPAR 

FILES 

FREQX 

GCOND 

HEACOM 

IONOSN 

MAGLNK 

OPTION 

ORIGIN 

PARAM 

PATH 

PHEN 

RANTEN 

SYSPAR 

TANTEN 

TIME 

USECB 

WEDEPO 

WOG 

WON 

WOX 

★ 

FXOPT  is  a Honeywell  system  subroutine  used  to  control  error 
messages  in  the  event  of  underflows  and  overflows. 


Hr* 

LENTRY  is  a Honeywell  system  subroutine  used  in  linking 
(overlaying). 


Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

ABFACR 

AIRY 

AIRY1 

FREQX 

GESCOM 

AIRY 

CABS 

AIRY1 

CEXP 

CLOG 

CSQRT 

CMPLXJ 

APPROX 

ABFACR 

CABS 

CMPLXJ 

FG 

CEXP 

CONST 

QSUBI 

CSQRT 

FREQX 

ROOT 

GESCOM 

SURFQ 

GND 

BHRHS 

ALOG 

C INPUT 

FREQX 

IONOSN 

PATH IN 

BUDDEN 

ENNU 

COS 

CONST 

ERINT 

SIN 

FIELD 

EMATRX 

FREQX 

ESTART 

INMM 

MMATRX 

RCCOM 

RINT 

SEBUG 

RSTART 

TCURVE 

CANG 

ATAN2 

CHEMDQ 

ATMOS U 

ALOG 

ATMOSN 

DRATE 

EXP 

ATMOST 

DTNEQ 

SQRT 

C INPUT 

IONOSU 

SPECDQ 

FDSRAT 

CONFAC 

COS 

CONST 

SIN 

GEOCOM 

SQRT 

PATH  IN 

TAN 

PROPO 

23 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

CONMAP 

COS 

CONST 

SIN 

DELTAT 

CMPLXJ 

FREQX 

GESCOM 

EHOP 

CONFAC 

CEXP 

CMPLXJ 

FMAG 

COVER1 

LENTRP 

COVER2 

REFLCO 

FREQX 

FVCOM 

GEOCOM 

OPTION 

OIJTREF 

SEBUG 

TCURVE 

EIGENV 

BUDDEN 

CABS 

CANGLE 

HELP 

CEXP 

CMPLXJ 

RGND 

CONST 

FREQX 

RGNDOU 

RPRIME 

SEBUG 

SYSTEM 

TCURVE 

EIGSRH 

EXCFAC 

CABS 

CANGLE 

HELP 

CCOS 

CMPLXJ 

HTGAIN 

CEXP 

CONST 

RGND 

CLOG 

DCOM 

WFDFDT 

COS 

EXCCOM 

WFZERO 

CSIN 

FREQX 

WINIFS 

CSQRT 

F1F2C 

W IN  I RB 

SQRT 

HCOM 

WSETRH 

»■ 

HTGMC 

INTEGR 

MODCOM 

OPTION 

PARAMC 

PATH 

RBCOM 

Table  2.  Continued. 


WEDCOM  IV  Routines 

C ' i 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

RGNDOU 

SIDESC 

SYSTEM 

TCURVE 

ELDEN2 

DATE 

LENTRY* 

ATMOUP 

JULIAN 

BTUBE 

LINKER 

CHMOVL 

RINTER 

CINPUT 

SOLORB 

CONST 

SOLZEN 

DBREG 

UN  I TV 

OPTION 

VECLIN 

ORIGIN 

VECM 

PROP 

XMAG 

ZTTOUT 

TIME 

ELFMOD 

BFIELD 

CEXP 

CANGLE 

BUDDEN 

CLOG 

CINPUT 

EIGENV 

CSQRT 

CMPLXJ 

GEOCOR 

CONST 

OUT  I ON 

ENEENP 

TEh 

FIELD 

FREQX 

GCOND 

GND 

10N0SN 

OPTION 

PATH 

RCCOM 

RGNDOU 

RPRIM 

SEBUG 

TANTEN 

TCURVE 

XMINS 

★ 

LENTRY  is  a Honeywell  system  subroutine  used  in  linking 
(overlaying) . 


Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

EMATRX 

CABS 

CMPLXJ 

EOUTM 

FIELD 

FREQX 

INMM 

SEBUG 

TCURVE 

XMINS 

ENNU 

ALOG 

CINPUT 

EXP 

ENEENP 

INMM 

ENVIRM 

ATMOSU 

ACOS 

ATMOS S 

AZELR 

COS 

BREG 

CONMAP 

LENTRY* 

CINPUT 

COORDV 

SIN 

COMINT 

COORDX 

SQRT 

CONST 

DATE 

DTUBE 

EDEPNB 

EDTOVL 

ELCOL 

FREQX 

ELDEN2 

GCOND 

GEOCOR 

IONOSN 

LINKER 

MAGLNK 

LOCLAX 

OPTION 

MAGFIT 

ORIGIN 

MAGFLD 

PATH 

MATMUL 

PRECAL 

SUBVEC 

PROP 

UNITV 

RANTEN 

VECM 

SOURCE 

XMAG 

TANTEN 

ZTTOUT 

TIME 

WRATE 

E PR I ME 

CEXP 

CANGLE 

EOUT 

FIELD 

FREQX 

RPRIM 

SEBUG 

★ 

LENTRY  is  a Honeywell  system  subroutine  used  in  linking  (over- 

laying). 

I 

Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

ERINT 

EMATRX 

CLOG 

CINPUT 

ENNU 

SQRT 

FIELD 

EPRIME 

OUTRST 

ESMAT 

RPRIM 

SEBUG 

TCURVE 

XMINS 

ESMAT 

CANGLE 

EOUT 

EOUTM 

FIELD 

SEBUG 

ESTART 

HELP 

CABS 

CANGLE 

CEXP 

CONST 

CLOG 

EOUTM 

CSQRT 

FIELD 

OUTRST 

SEBUG 

OTGHN 

CEXP 

CMPLXD 

FREQX 

HGCOM 

EXCFAC 

CSQRT 

CANGLE 

EXCCOM 

RGNDOU 

FEDELA 

AS  IN 

CONST 

COS 

PROPO 

SIN 

FINDHP 

FEDELH 

ALOG 

PATHIN 

COS 

FINGEO 

FDELH 

COS 

CONST 

FINDHP 

SIN 

FGNDC 

TAN 

GEOCOM 

PATHIN 

PROPO 

Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

FG 

AIRY 

AT  AN  2 

AIRY1 

CABS 

CMPLXJ 

CEXP 

CLOG 

CSQRT 

CONST 

FMAG 

AIRY 

CABS 

AIRY1 

GREFL 

CEXP 

CMPLXJ 

RUNG 

COS 

CONST 

SURFQ 

SQRT 

COVER1 

FGNDC 

FREQV 

FVCOM 

GEOCOM 

PATH  IN 

PROPO 

GNDWV 

FMAG 

FVCOM 

GCOND 

OPTION 

PATH IN 

GREFL 

COS 

CSQRT 

SIN 

PROPO 

GUESS 

ABFACR 

CLOG 

AIRY1 

AIRY 

CSQRT 

CANGLE 

APPROX 

CMPLXJ 

ISOREF 

CONST 

QSUBI 

EXCCOM 

SURFQ 

FREQX 

GESCOM 

GND 

MODCOM 

RGNDOU 

SYSTEM 

TCURVE 

Table  2.  Continued 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

HTGAIN 

AIRY 

CEXP 

AIRY1 

CSQRT 

CANGLE 

EXP 

CMPLXJ 

CONST 

EXCCOM 

FREQX 

GND 

HTGMC 

RGNDOU 

TCURVE 

IONOSU 

RATE 

EXP 

SIN 

SQRT 

INPUT 

AZELR 

EXIT 

BURST 

COORDX 

CINPUT 

GEOCOR 

CONST 

HELP 

DEVICE 

LOCLAX 

FILES 

MATMUL 

FREQX 

UANTD 

GCOND 

VECLIN 

OPTION 

VECM 

ORIGIN 

Xl'iAG 

PATH 

WOG1 

RANTEN 

WON1 

SYSPAR 

WOX1 

TANTEN 

USECB 

ISOREF 

CANG 

ALOG 

CINPUT 

ENNU 

CABS 

CMPLXJ 

CEXP 

CONST 

CSQRT 

FREQX 

EXP 

GESCOM 

INMM 

SEBUG 

TCURVE 

29 


Table  2.  Continued. 


Routine 

WEDCOM  IV  Routines 
Called  Directly 

Library  Routines 
Called  Directly 

Common  Blocks 
Used 

ISORLF 

LENTRP 

ENNU 

ALOG 

CEXP 

CSQRT 

EXP 

CINPUT 

CMPLXJ 

FREQX 

GESX 

INMM 

TCURVE 

CINPUT 

CONST 

ENEENP 

GCOND 

IONOSN 

PATH  IN 

LFHOP 

AZELR 

ALOGIO 

CINPUT 

BFIELD 

CABS 

CONST 

BHRHS 

CANG 

COVER1 

COORDX 

COS 

COVER3 

EHOP 

CSQRT 

ENEENP 

FINGEO 

EXP 

FIELD 

GEOCOR 

SIN 

FREQX 

GNDWV 

OUT  I ON 

TFIELD 

UNITV 

VECFi 

XMAG 

SQRT 

FVCOM 

GCOND 

IONOSN 

OPTION 

PATH 

PATHIN 

PROPO 

SEBUG 

SYSPAR 

TANTEN 

VHFLD 

XM  INS 

I'lATINV 

MDHNKL 

HELP 

CABS 

SQRT 

CABS 

CEXP 

CSQRT 

Table  2. 
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\ 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

MATRX 

CABS 

CMPLXJ 

FIELD 

FREQX 

INMM 

OUTMM 

TCURVE 

XMINS 

NODCON 

HTGAIN 

AT  AN  2 

ACOM 

MATIN V 

CABS 

AIRY1 

CANG 

CANGLE 

CCOS 

CMPLXJ 

CEXP 

CONST 

CSIN 

FREQX 

HTGMC 

MODCOM 

OPTION 

TCURVE 

VLFCOM 

MODEZ 

EIGENV 

CABS 

CANGLE 

EXCFAC 

CEXP 

CMPLXJ 

GUESS 

CONST 

HTGAIN 

EXCCOM 

FREQX 

HTGMC 

MODCOM 

OPTION 

PATH 

RGNDOU 

SYSTEM 

TCURVE 

OUT  I ON 

ENNU 

CSQRT 

CINPUT 

CMPLXJ 

CONST 

FIELD 

FREQX 

INMM 

OPTION 

TCURVE 

XMINS 

31 


Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

QSUBI 

CEXP 

CMPLXJ 

CSQRT 

CONST 

FREQX 

GESCOM 

REFLCO 

ENNU 

CEXP 

CANGLE 

ISORLF 

COS 

CINPUT 

MMATRX 

CSQRT 

CMPLXJ 

RINT 

SIN 

CONST 

RSTART 

FIELD 

FREQX 

GESX 

INMM 

OPTION 

OUTREF 

PARAM 

PROPO 

RPRIM 

TCURVE 

XMINS 

RGND 

AIRY 

CSQRT 

AIRY1 

CANGLE 

CONST 

FREQX 

GND 

RGNDOU 

SEBUG 

TCURVE 

RINT 

ENNU 

SQRT 

CINPUT 

MATRX 

FIELD 

RPRIME 

OUTRST 

SMAT 

RPRIME 

TCURVE 

XMINS 

ROOT 

ABFACR 

DELTAT 

HELP 

QSUBI 

CABS 

Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

RPR I ME 

FIELD 

FREQX 

OUTS 

RPRIM 

RSTART 

CABS 

CANGLE 

CEXP 

CONST 

CLOG 

FIELD 

CSQRT 

OUTMM 

OUTRST 

RUNG 

SMAT 

CANGLE 

FIELD 

OUTMM 

OUTS 

SURQ 

CSQRT 

CMPLXJ 

FREQX 

TEM 

AZELR 

CABS 

CINPUT 

COORDX 

CANG 

CMPLXJ 

ETGAIN 

CEXP 

CONST 

GEOCOR 

COS 

FREQX 

CSQRT 

HGCOM 

SIN 

OPTION 

SQRT 

PATH 

SYSPAR 

TANTEN 

TFIELD 

CABS 

CMPLXJ 

CEXP 

C0VER2 

CLOG 

COVER3 

SQRT 

GEOCOM 

OPTION 

VHFLD 

UANTD 

CINPUT 

TANTEN 

33 


Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

VLFMCV 

MODCON 

ACOS 

ACOM 

AL0G10 

CINPUT 

ATAN 

CMPLXJ 

CABS 

CONST 

CANG 

FREQX 

CEXP 

MODCOM 

CSIN 

OPTION 

EXP 

PATH 

SIN 

SYSPAR 

SQRT 

SYSTEM 

VLFCOMl 

VLFMOD 

AZELR 

COS 

ACOM 

BFIELD 

LENTRY* 

CINPUT 

COORDX 

SIN 

CONST 

GEOCOR 

1 

ENEENP 

LINKER 

FIELD 

OUT  I ON 

FREQX 

XMAG 

GCOND 

GND 

HCOM 

IONOSN 

MODCOM 

OPTION 

PARAM 

PATH 

RCCOM 

SEBUG 

SYSPAR 

SYSTEM 

TANTEN 

TCURVE 

XMINS 

VLFCOM 

* 


s a 


s 


s 


ie  use* 
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Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Routine 

Called  Directly 

Called  Directly 

VLFWKB 

ACOS 

ALOGIO 

ATAN 

CABS 

CANG 

CEXP 

CSIN 

CSQRT 

EXP 

SIN 

SQRT 

WDNATM 

ATMOSU 

ALOG 

CHEMDQ 

IONOSU 

RATE 

SQRT 

WFCTVL 

WLGRNG 

WRBARS 

WFDFDT 

WLGDER 

WRBARS 

CSIN 

WFINDF 

WFCTVL 

WFSINT 

MDHNKL 

ALOG 

CCOS 

CEXP 

EXP 

Common  Blocks 
Used 


CINPUT 

CMPLXJ 

CONST 

FREQX 

MODLOM 

OPTION 

PATH 

SYSPAR 

SYSTEM 

VLFCOM 


ATMOSN 

ATMOSS 

ATMOST 

ATMOUP 

CINPUT 

IONOSN 

OPTION 

OUT 

IDERVC 

INTEGR 

RBCOM 

THETAC 

CONST 

IDERVC 

INTEGR 

RBCOM 

THETAC 

TMCCOM 

CANGLE 

CMPLXJ 

CONST 

FREQX 

HCOM 

IDERVC 

INTEGR 

THETAC 


Table  2.  Continued. 


WEDCOM  IV  Routines 

Library  Routines 

Common  Blocks 

Routine 

Called  Directly 

Called  Directly 

Used 

WFZERO 

HELP 

TMCCOM 

WFINDF 

WNOMES 

WQUAD 

ZLSLOM 

WLGRNG 

CCOS 

CONST 

INTEGR 

THETAC 

WLGINT 

WNOMES 

HELP 

WFDFDT 

WFINDF 

WQUAD 

CABS 

ZLSCOM 

WQUAD 

SQRT 

WRBARS 

MDHNKL 

ALOG 

CMPLXJ 

CCOS 

CONST 

CEXP 

DCOM 

CSQRT 

F1F2C 

EXP 

FREQX 

GND 

HCOM 

IDERVC 

RBCOM 

THETAC 

WSETRH 

BUDDEN 

CABS 

CANGLE 

HELP 

CCOS 

CONST 

WFSINT 

CSIN 

DCOM 

WINIFS 

IDERVC 

WINILG 

INTEGR 

WLGDER 

PARAMC 

WSETLG 

RPRIM 

XFER 

SIDESC 

TCURVE 

THETAC 

WLGINT 

XFER 


SECTION  2 

ROUTINE  DESCRIPTIONS 


CONTROL  ROUTINE 

This  routine  is  the  main  driver  routine  for  WEDCOM  IV.  A 
simplified  flowchart  for  the  routine  is  shown  in  Figure  1. 

Before  starting  the  problem,  parameters,  constants,  and  in- 
put default  values  are  defined.  The  problem  loop  refers  to  a given 
input  specification  and  the  problem  number  is  changed  each  time  con- 
trol is  transferred  to  the  input  module  and  new  input  data  obtained. 

After  input  data  are  obtained  a check  is  made  to  see  if  phe- 
nomenology or  environmental  calculations  have  been  previously  made 
for  the  specified  input  conditions.  If  phenomenology  or  environ- 
mental calculations  are  required,  a check  is  made  to  see  if  the  ambi- 
ent atmosphere  and  magnetic  dipole  field  parameters  are  to  be  calcu- 
lated at  only  the  input  origin  (input  option).  If  they  are  (MAMB  = 0) 
the  ambient  atmosphere  and  ambient  magnetic  field  routines  are  ini- 
tialized and  the  time  of  day  at  the  input  origin  determined. 

Next,  a loop  over  calculation  times  is  started  and  the  phenome- 
nology module  (driver  routine  PHENOM  for  the  RANC  IV  phenomenology  or 
driver  routine  PHENOR  for  the  ROSCOE  phenomenology)  is  called  to  ob- 
tain fireball  and  debris  quantities.  If  the  atmospheric  wind  model 
is  to  be  used  (input  option),  the  wind  module  (driver  routine  DBSTAM) 
is  called  to  determine  wind  effects  on  debris  regions.  If  the  ROSCOE 
phenomenology  is  used,  the  low-altitude  debris  regions  are  converted 
to  equivalent  RANC  debris  regions.  This  is  done  in  the  current  ver- 
sion of  the  WEDCOM  code  to  simplify  geometry  calculations  required 
when  determining  ionization  and  absorption.  Use  of  the  ROSCOE 


low-altitude  model  retains  multiburst  effects  and  simplification  of 
the  detailed  geometry  (eg,  modified  Oval  of  Cassini)  is  not  critical. 
To  reduce  storage  requirements  the  debris  region  data  are  stored  in 
a weapon  file.  Nominal  and  detailed  outputs  describing  the  debris 
regions  are  also  written  out. 

After  the  debris  region  data  have  been  determined  environ- 
mental quantities  between  the  transmitter  and  receiver  are  found  by 
calling  the  environment  module  (driver  routine  ENVIRM) . Calculations 
in  the  environment  module  are  made  at  specified  altitudes  on  vertical 
paths  uniformly  spaced  along  the  great  circle  path  between  transmitter 
and  receiver.  Calculations  are  made  for  all  calculation  times  and 
the  results  stored  in  a file  to  reduce  storage  requirements.  The 
environmental  data  are  then  reordered  so  that  the  data  are  available 
for  all  paths  at  each  calculation  time  (order  necessary  to  determine 
propagation  effects)  by  calling  routine  REORDR. 

Propagation  effects  are  found  by  starting  a frequency  loop 
and  calling  the  ELF,  VLF,  or  LF  propagation  modules.  For  ELF  and 
VLF,  calculations  are  made  within  the  module  for  either  a single 
receiver  location  or  for  multiple  receiver  locations,  depending  on 
user  input.  In  the  LF  propagation  module,  calculations  are  made  for 
a single  receiver  location.  A loop  over  multiple  receiver  locations 
is  made  in  the  control  routine  if  multiple  receiver  location  data 
are  requested. 
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SUBROUTINE  ABFACR 

This  routine  calculates  terms  in  the  isotropic  VLF  mode 
equat  ion . 

' 

The  terms  A and  B,  defined  in  Equations  3-6  and  3-7  in 
Section  3,  Volume  3,  are  computed.  In  addition,  the  calculated 
values  for  some  parameters  are  saved  for  later  use  in  routine 


DELTAT. 


SUBROUTINE  AIRY 


I 


f 
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This  routine  calculates  Airy  function  values  and  their 
derivatives . 

The  Airy  functions  as  defined  by  Wait  (Reference  3)  can  be 

written 


Wx(t) 

= ay^t)  + by2(t) 

(1) 

w2(T) 

= a*y1(t)  + b*y2(t) 

(2) 

Wj(  t) 

= ay{(t)  + by 2 (t) 

(3) 

w'(t) 

= a*y'(t)  + b*y2(t)  , 

(4) 

where 

a = 1 

.089929069  - i 0.629270841 

(5) 

b = 0 

.794570425  + i 0.458745449 

(6) 

and  the  asterisk  indicates  complex  conjugate. 

The  functions  y ^ (t ) 

and  >'-,(t)  are 

infinite  series  in  the  argument 

(t) , and  the  prime 

indicates  derivative  with  respect  to  the  argument.  The  y functions 

can  be  written  in  a form  suitable  for  computation  as  follows: 

yl  = 

in 

1 + L yln 
n=l 

(7) 

t'J 

II 

■(■ ' &•) 

(8) 

t 

yl  = 

•H 

(9) 

v. 

y 2 = 

m 

1 *!>(„ 
n=l 

(10) 
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where 


^0  y20  yl0  y20  1 


(11) 


and  the  additional  terms  are  obtained  from  the  recursion  relations 
yln 


yl,n-l 

1 

Pt 

0-> 

(12) 

3n(3n  - 1) 

— \r 

1 

t3 

(13) 

y2,n-l 

3n(3n  + 1) 

= v' 

1 

t3 

(14) 

yl,n-l 

3n(3n  + 2) 

= y2 ,n- 1 

1 

t3 

(15) 

3n(3n  - 2) 

number  of  terms,  m. 

required 

for  a given  value  of  the 

2n 


argument,  t,  is  defined  so  that  the  magnitude  of  the  mth  term  in  any 

Q 

of  the  series  7 through  10  is  less  than  10~  . The  number  of  terms 
required  for  the  desired  accuracy  as  a function  of  the  magnitude  of 
t can  be  approximated  by  a straight  line  defined  by 


m = 4 + 2,8  t 


(16) 


I 1 
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SUBROUTINE  APPROX 


This  routine  calculates  initial  values  for  the  approximate 
mode  equation.  A flowchart  for  the  routine  is  shown  in  Figure  2. 

Solutions  to  the  approximate  mode  equation  are  obtained  em- 
ploying the  Newton  iteration  algorithm  after  an  appropriate  guess  is 
calculated.  Four  approximate  solutions  are  used  in  obtaining  the 
guess  (see  Section  3,  Volume  3).  The  solution  starts  with  a modified 
flat-earth  approximation  (Branch  1)  for  the  mode  angle.  If  the  flat- 
earth  approximation  cannot  be  used,  a test  is  made  to  see  if  it  can 
be  used  as  a first  guess  for  the  iterative  solution  of  the  mode  equa- 
tion. If  this  test  fails,  a grazing-earth  incident-angle  approxima- 
tion (Branch  2)  is  computed  and  tested.  The  first  two  branches  pro- 
vide an  adequate  guess  for  modes  greater  than  1.  If  neither  of  these 
approximations  is  good  for  a mode  1 solution,  a third  guess  is  ob- 
tained assuming  that  the  mode  belongs  to  the  "whispering  gallery" 
type  (Branch  3).  If  this  third  calculation  does  not  provide  a good 
guess,  either  the  best  guess  from  the  first  three  branches  is  found 
or  a search  through  a matrix  of  initial  values  is  made  (Branch  4), 
where  the  matrix  elements  belong  to  the  transition  region  between 
whispering  gallery  and  grazing  incidence  approximations.  The  Branch 
4 calculation  is  used  only  for  mode  1 vertical  field  calculations, 
since  the  first  three  approximations  normally  provide  a good  guess 
for  higher  order  vertical  or  horizontal  field  calculations. 

After  the  best  guess  is  determined,  the  solution  of  the  mode 
equations  is  calculated  employing  routine  ROOT.  Occasionally,  for 
the  smaller  ground  conductivities,  the  iteration  technique  in  routine 
ROOT  does  not  converge  to  a proper  solution.  If  a good  solution  is 
not  obtained,  an  increased  ground  conductivity  is  chosen  for  which  a 
mode  solution  can  be  accurately  determined.  The  code  then  returns 
to  the  flat-earth  approximation,  and  the  process  described  above  is 
repeated  to  determine  the  best  guess  for  this  new  conductivity.  The 
solution  for  the  mode  equation  is  determined  (Branch  5)  and  used  as 


f 


a guess  for  a reduced  conductivity.  This  process  is  repeated  until 
the  conductivity  is  equivalent  to  the  original  specified  value. 


SUBROUTINE  BUDDEN 


This  routine  is  the  driver  routine  for  ELF  and  VLF  reflection 
coefficient  calculations  in  an  anisotropic  ionosphere. 

The  reflection  coefficients  are  found  by  dividing  the  iono- 
sphere profile  into  slabs  and  then  performing  a numerical  integration 
of  differential  equations  for  the  reflection  coefficients  starting 
at  the  top  slab. 

First,  magnetic  field  terms  are  determined  and  routine  ENNU 
called  to  evaluate  the  ionization  and  collision  frequencies  at  the 
top  slab.  Then  the  starting  values  for  the  reflection  coefficients 
are  found  by  calling  routines  MMATRX  and  RSTART  (VLF  propagation)  or 
routines  EMATRX  and  ESTART  (ELF  propagation).  After  the  starting 
values  are  determined  the  numerical  integration  is  performed  by 
calling  routine  RINT  (VLF  propagation)  or  ERINT  (ELF  propagation). 


FUNCTION  CANG 


This  routine  calculates  the  argument  of  a complex  number  (in 


radians) 


The  argument  is  computed  from: 

cang  = tan_1  Iris  radians  > 

where  Re  and  Im  are  the  real  and  imaginary  parts,  and  ARC  is  the 
complex  number. 
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SUBROUTINE  CHEMDQ 

This  routine  is  a driver  routine  to  determine  steady-state 
ionization  and  electron  and  ion  collision  frequencies  for  a given 
delayed  ionization  source  or  an  ambient  ionization  source  at  a speci- 
fied D-region  altitude.  The  routine  is  the  same  as  routine  CHEMDQ 
used  in  WEP1I  VI  except  for  the  ion-neutral  collision  frequency 
calculation. 

First,  routine  ATMOSU  is  called  to  obtain  the  ambient  neutral 
species  at  the  altitude  of  interest  and  routine  DRATE  is  called  to 
determine  the  values  of  the  reaction  rate  coefficients.  Then,  if 
the  input  value  of  the  ion-pair  production  rate  (Q)  is  zero,  routine 
IONOSU  is  called  to  obtain  the  ambient  ion-pair  production  rate. 

Next,  routine  SPECDQ  is  called  to  determine  the  change  in 
neutral  species  and  routine  DTNEQ  is  called  to  determine  the  electron 
and  ion  densities  due  to  the  ion-pair  production  rate.  The  final 
calculations  are  for  the  electron  and  ion  collision  frequencies.  A 
fit  to  an  expression  developed  at  the  Naval  Ocean  Systems  Center 
(Reference  4)  is  used  for  the  ion-neutral  collision  frequency. 


SUBROUTINE  CONFAC 


■ 

! 

This  routine  computes  the  ionospheric  convergence  factor. 

The  ionospheric  convergence  factor  is  computed  for  either  a 
uniform  or  nonuniform  ionosphere  along  the  propagation  path.  For 
either  situation,  applicable  approximate  representations  are  used 
for  the  convergence  factor  when  geometric  optics  apply  or  when  dif- 
fractive corrections  are  required.  The  convergence  factor  is  com- 
puted by  evaluating  equations  given  in  Section  4,  Volume  3. 


SUBROUTINE  CONMAP 


l! 


This  routine  provides  ground  conductivity  for  a specified 
geographic  location.  It  contains  a modified  version  of  the  Institute 
of  Telecommunication  Sciences  world  conductivity  map. 

Fourier  analysis  is  applied  to  the  longitudinal  variation  in 
conductivity,  while  the  latitudinal  variation  is  represented  by  a 
weighted  least  squares  polynomial  in  the  sine  of  the  latitude. 

After  the  initialization  of  parameters,  the  latitudinal 
values  of  the  Fourier  coefficients  are  determined.  The  Fourier  coef- 
ficients are  used  with  the  given  longitude  value  to  determine  a value 
that  is  directly  related  to  the  conductivity  of  the  desired  location. 


SUBROUTINE  EHOP 


This  routine  computes  vertical  and  horizontal  field  strengths 
for  individual  skywaves.  A flowchart  for  the  routine  is  shown  in 
Figure  4. 

Routine  EHOP  is  called  after  all  normal  and  disturbed  iono- 
spheric parameters  and  the  ray  path  geometry  have  been  established. 

Two  major  calculation  loops  are  used.  In  one,  field  strengths  are 
obtained  for  skywaves  with  a geometry  that  requires  diffractive  cor- 
rections. In  the  other,  field  strengths  are  calculated  for  skywaves 
where  geometric  optics  approximations  are  appropriate. 

In  both  loops  interpolation  procedures  are  used  to  define 
ionospheric  and  ground  parameters  at  the  ionospheric  or  ground  reflec- 
tion points.  In  the  first  loop,  diffractive  correction  terms  in  the 
vicinity  of  ground  reflection  points  are  defined. 

An  inner  loop  is  also  used  for  geometric  optics  approximations 
when  either  transmitter  or  receiver  antennas  are  elevated.  In  this 
case,  fields  for  two  distinct  rays  with  an  equ-’l  number  of  iono- 
spheric reflections  are  computed  and  summed. 

The  actual  equations  evaluated  and  the  procedures  used  for 
defining  the  number  of  diffractive  paths  are  discussed  in  Section  4, 
Volume  3. 
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Figure  4.  Flowchart  for  subroutine  EHOP. 


SUBROUTINE  EIGENV 
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This  routine  calculates  reflection  coefficients  and  the  eigen- 
value and  derivative  of  the  mode  equation  using  the  Newton-Raphson 
procedure.  A flowchart  for  the  routine  is  shown  in  Figure  5. 

The  initialization  of  iteration  parameters  is  followed  by  a 
second-order  1)0  loop  that  determines  an  improved  eigenvalue  for  the 
initial  eigenvalue  guess  (an  initial  guess  given  as  input  is  used  as 
the  starting  value)  and  the  incremented  value  of  the  eigenvalue. 

Both  of  these  calculations  are  made  for  the  first  iteration.  For 
subsequent  iterations  one  or  both  are  done,  depending  on  the  computed 
eigenvalue  increment.  An  input  option  determines  whether  the  value 
of  the  exact  mode  equation  is  to  be  computed,  or  if  a vertical  or 
horizontal  polarization  approximation  is  to  be  made.  If  the  process 
diverges  greatly  from  the  initial  guess,  the  eigenvalue  is  set  to  the 
initial  guess  and  control  returns  to  the  calling  routine. 

The  value  of  the  derivative  to  the  mode  equation  is  found, 
and  a test  is  made  to  determine  if  the  eigenvalue  has  been  determined. 
If  so,  control  returns  to  the  calling  routine.  Otherwise,  the  eigen- 
value increment  for  the  next  iteration  is  calculated,  and  the  itera- 
tion index  is  updated.  Control  returns  to  the  calling  routine  if 
there  are  more  than  20  iterations.  If  not,  tests  are  made  to  deter- 
mine if  the  convergence  criterion  has  been  satisfied.  If  these  tests 
have  been  satisfied,  a test  is  made  to  see  if  the  proper  value  of  the 
derivative  of  the  mode  equation  was  calculated.  If  not,  it  is  calcu- 
lated before  returning  control  to  the  calling  routine. 

If  another  iteration  is  necessary,  an  eigenvalue  increment  is 
calculated.  If  this  increment  is  not  larger  than  a predetermined 
value,  one  of  the  prior  mode  equation  values  is  used  and  only  one  new 
calculation  of  the  mode  equation  is  made.  Otherwise,  two  new  values 
of  the  mode  equation  are  computed  in  the  determination  of  the 
derivat ive. 


The  eigenvalue  and  the  derivative  of  the  mode  equation  are 
referenced  to  altitude  l).  Altitude  D is  set  to  0 km  for  ELI-'  calcula- 
tions and  for  VLF  calculations  that  do  not  use  the  Naval  Ocean  Systems 
Center  mode  search  models  (see  description  of  routine  EIGSRH) . The 
reflection  coefficients  are  referenced  to  altitude  HBOT.  Altitude 

HBOT  may  be  higher  than  altitude  D for  VLF  propagation.  * 


SUBROUTINE  EIGSRH 
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This  routine  computes  solutions  to  the  mode  equation  for  a 
given  ionosphere  profile  and  ground  conductivity  from  a mode  search 
routine  developed  by  the  Naval  Ocean  Systems  Center  (Reference  5). 
Height  gain  and  excitation  factors  are  also  calculated.  A simplified 
flowchart  for  the  routine  is  shown  in  Figure  6. 

Solutions  to  the  mode  equation  (eigenangles)  are  found  using 
Naval  Ocean  Systems  Center  mode  search  routines  to  search  for  solu- 
tions within  a specified  region  (rectangle)  in  the  angle  of  incidence 
complex  plane.  A total  of  three  rectangles  are  defined  and  a maxi- 
mum of  15  eigenangles  are  found.  Of  these,  only  the  10  most  impor- 
tant are  saved  for  use  in  determining  received  field  strengths. 

For  each  rectangle  in  the  angle  of  incidence  complex  plane 
routine  WHZERO  is  called  to  determine  eigenangles.  The  eigenangles 
are  ordered  in  terms  of  their  real  part  and  solutions  outside  the 
current  rectangle  are  discarded.  Routines  EXCFAC  and  HTGAIN  are 
called  to  determine  the  excitation  and  height  gain  factors. 

Next,  a weighted  attenuation  rate  that  includes  the  effect  of 
excitation  factors  at  the  transmitter  and  for  the  ionospheric  profile 
for  the  current  vertical  path  is  computed.  The  attenuation  rate  is 
used  to  select  up  to  10  modes.  After  the  modes  have  been  selected 
they  are  reordered  so  that  the  real  part  of  the  eigenangles  decrease 
as  the  mode  number  increases.  Output  is  prepared  describing  the 
modes  and  if  detailed  output  is  requested  (input  option),  details  of 
the  mode  characteristics  height  gain  factors,  and  excitation  factors 
are  written  out. 


* 
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ANGLE  OF  INCIDENCE  COMPLEX 
PLANE  REGION  LOOP 


Figure  6.  Flowchart  for  subroutine  EIGSRH. 
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SUBROUTINE  ELCOL 

This  routine  computes  the  colatitude  and  east  longitude  of 
point  2,  given  the  colatitude  and  east  longitude  of  point  1 and  the 
azimuth  (radians)  and  ground  range  (km)  of  point  2 from  point  1.  The 
supplement  of  the  azimuth  of  point  1 from  point  2 is  also  found. 

The  azimuth  of  point  2 from  point  1 is  positive,  measured 
clockwise  from  north.  The  supplement  of  the  azimuth  of  point  1 from 
point  2 is  determined  as  a positive  number. 

The  colatitude  of  point  2 is  obtained  from 

cos  $ , = cos  4^  cos  R/Re  + sin  sin  R/R^  cos  a , (1) 


where 


<4> ^ = colatitude  of  point  2 
<J>j  = colatitude  of  point  1 

R = ground  range  between  point  1 and  point  2 
e = earth's  radius 


A12 

A12s  n 

A 

C\l 

•-H 

< 

= azimuth  of  point  2 with  respect  to  point  1. 
The  east  longitude  of  point  2 is  obtained  from 


0J  + A0 


A12  5 * 


0J  + 2 it  -A0  Aj2  > IT 
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The  supplement  of  the  azimuth  of  point  1 with  respect  to  point  2 i; 
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SUBROUTINE  ELDEN2 


This  routine  is  the  driver  routine  for  ionization  and  colli- 
sion frequency  calculations  at  a specified  location.  The  routine  is 
the  same  as  routine  ELDEN2  used  in  the  WEPH  VI  code  except  for  the 
collision  frequency  calculations. 

For  altitudes  at  or  below  100  km  routine  CHEMD2  is  called  to 
obtain  ionization  and  collision  frequency  quantities.  For  altitudes 
above  100  km  a calculation  time  loop  is  set  up  and  routine  ENEF  called 
for  each  calculation  time.  Routine  LINKER  is  called  to  overlay  the 
D-region  (CHEMD2)  and  E-  and  F-region  (ENEF)  routines.  The  labeled 
common  block  CHMOVL  is  used  to  transfer  input  to  and  obtain  output 
from  routines  CHEMD2  and  ENEF. 

For  altitudes  above  100  km  a calculation  time  loop  is  set  up 
and  for  each  calculation  time  the  time  of  day  (day  or  night)  is  com- 
puted. For  the  first  calculation  altitude  intersection  altitudes  of 
the  ray  path  with  the  beta  tube  are  computed  for  each  debris  region. 
The  intersection  altitudes  are  stored  on  a data  file  for  use  at  sub- 
sequent calculation  altitudes.  If  detailed  ionization  output  is  re- 
quested (input  option),  data  is  obtained  from  file  LFIDO  and  written 
out  on  file  LFDO. 


SUBROUTINE  ELFMOD 


This  routine  is  the  ELF  propagation  driver  routine  and  de- 
termines anisotropic  reflection  coefficients,  mode  solutions,  and 
excitation  factors.  Only  the  transverse  electromagnetic  wave  is 
analyzed.  A simplified  flowchart  for  the  routine  is  shown  in  Figures 
7 and  8. 

First  reflection  coefficients,  excitation  factors,  and  height 
gain  parameters  are  initialized.  Then  a time  loop  over  the  calcula- 
tion times  is  started  and  an  ambient  calculation  flag  (determines 
whether  calculations  will  be  made  for  ambient  conditions)  is  set  from 
input  parameters. 

Next,  a loop  over  the  vertical  paths  between  transmitter  and 
receiver  is  started.  For  each  path,  magnetic  field  and  ground  con- 
ductivity quantities  are  established  and  the  level  of  detailed  output 
requested  determined.  The  ionization  and  collision  frequencies  along 
the  vertical  path  for  disturbed  conditions  are  obtained  from  the 
environment  file  and  routine  OUTION  called  to  determine  the  top  and 
bottom  altitudes  to  be  used  in  computing  reflection  coefficients. 

For  ELF  propagation  there  are  NP-2  vertical  paths  defined 
between  the  transmitter  and  receiver;  NPT2  paths  along  the  short 
great  circle  path  and  the  rest  along  the  long  great  circle  path  (only 
defined  if  calculations  along  the  long  great  circle  path  are  to  be 
made).  The  last  two  paths  (L  = LP  - 1 and  L = LP)  are  located  a 
Fresnel  zone  distance  normal  to  the  short  great  circle  path  between 
transmitter  and  receiver  and  are  used  to  determine  whether  the  ioni- 
zation normal  to  the  great  circle  propagation  path  is  spherically 
stratified  (an  assumption  used  in  the  propagation  calculations).  If 
the  attenuation  rate  (dB  per  1000  km)  for  either  of  the  last  paths 
is  less  than  one  half  the  attenuation  rate  at  the  path  midpoint,  a 
diffraction  flag  is  set  to  indicate  that  the  propagation  calculations 
may  be  in  error. 
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Calculations  for  the  reflection  coefficient  and  mode  solu- 
tions are  made  for  vertical  paths  between  the  transmitter  and  receiver 
(L  < NP  - 2) . If  the  input  option  to  minimize  the  number  of  reflec- 
tion coefficient  calculations  is  exercised  (MRCOEF  = 1),  tests  are 
first  made  to  see  how  much  the  top  altitude  obtained  from  routine 
OUTION,  the  magnetic  field  properties,  and  the  ground  conductivity 
have  changed  from  the  values  computed  for  the  previous  vertical  path. 

If  they  have  not  changed  by  more  than  specified  amounts,  the  eigen- 
angle,  ground  reflection  coefficient,  and  excitation  factors  previ- 
ously determined  are  used  for  the  current  path. 

Otherwise,  routine  BUDDEN  is  called  and  calculations  made  to 
obtain  initial  guesses  for  the  reflection  coefficients  and  eigenangle 
and  routine  EIGENV  called  to  obtain  final  values  (see  Figure  8).  If 
the  multiple  receiver  location  option  is  exercised  (KALCMR  = 1) , ground 
reflection  coefficients  and  excitation  factors  are  found  for  each  ver- 
tical path  on  the  short  great  circle  path  between  transmitter  and  re- 
ceiver. Otherwise,  ground  reflection  coefficients  and  excitation 
factors  are  only  found  for  the  vertical  paths  above  the  transmitter 
and  receiver  terminals.  If  field  strength  calculations  for  propaga- 
tion along  both  the  short  and  long  great  circle  paths  are  to  be  made, 
ground  reflection  coefficients  and  excitation  factors  at  the  trans- 
mitter and  receiver  terminals  are  computed  for  both  propagation  direc- 
tions. The  ground  reflection  coefficients  and  excitation  factors  for 
the  long  great  circle  path  direction  are  not  computed  for  receiver 
locations  between  the  transmitter  and  receiver  terminals  specified 
in  input. 

Detailed  output  describing  the  eigenangle,  attenuation  rate, 
relative  phase  velocity,  ionosphere  reflection  coefficients,  ground 
reflection  coefficients,  and  excitation  factors  are  written  out  on 
the  detailed  output  file  if  requested  (input  option). 
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SUBROUTINE  EMATRX 


This  routine  calculates  elements  of  the  ionospheric  suscepti- 
bility matrix,  including  magnetic  field  effects  and  earth  curvature. 

The  squares  of  the  ratios  of  electron,  positive  ion,  and  nega- 
tive ion  plasma  frequencies  to  the  signal  frequency  are  defined.  Next, 
the  corresponding  ratios  of  collision  frequencies  to  the  signal  fre- 
quency are  defined.  The  susceptibility  matrix  is  initialized. 

The  susceptibility  matrix  is  calculated  by  separately  calcu- 
lating and  then  summing  the  contributions  due  to  electrons  and  posi- 
tive and  negative  ions.  To  reduce  numerical  problems,  the  matrix 
quantities  are  normalized  by  the  square  of  the  ratio  of  electron 
plasma  frequency  to  signal  frequency.  If  magnetic  field  effects  are 
not  important,  only  the  diagonal  terms  are  calculated.  The  effects 
due  to  earth  curvature  are  included  by  an  artificial  modification  of 
the  diagonal  terms  in  the  susceptibility  matrix. 


SUBROUTINE  ENNU 


This  routine  determines  electron  and  ion  densities  and  colli- 
sion frequencies  at  specific  altitudes. 

Electron  and  ion  densities  and  collision  frequencies  are 
determined  by  exponential  interpolation  between  field  points.  The 
field  points  and  ionospheric  parameters  are  predetermined  and  avail- 
able through  labeled  common  blocks  CINPUT  and  ENEENP. 


SUBROUTINE  ENVIRM 


This  routine  is  the  driver  routine  for  propagation  environ- 
ment calculations.  A simplified  flowchart  for  the  routine  is  shown 
in  Figure  9. 

First,  the  number  of  vertical  paths  to  be  located  on  the  short 
great  circle  path  between  the  transmitter  and  receiver  terminals  are 
determined  and  the  location  of  the  vertical  path  closest  to  the  mid- 
point of  the  great  circle  path  found.  For  ELF  propagation  two  addi- 
tional paths  located  a Fresnel  zone  distance  normal  to  the  great 
circle  path  midpoint  are  defined.  These  paths  are  for  use  in  deter- 
mining a diffraction  flag.  Then,  the  ambient  atmosphere  and  magnetic 
routines  are  initialized  at  the  short  great  circle  path  midpoint. 

If  calculations  are  to  be  made  for  the  long  great  circle  path  be- 
tween transmitter  and  receiver  (KALCLP  > 0);  the  number  of  vertical 
paths  to  be  used  along  the  long  great  circle  path  are  determined. 

Next,  the  ambient  atmospheric  properties  at  the  short  great  circle 
path  midpoint  are  found  by  calling  routine  WDNATM. 

The  path  loop  over  the  vertical  paths  to  be  used  in  determi- 
ning the  environment  is  started  and  the  path  geometry  found.  If  the 
environment  above  100  km  is  to  be  determined  (ELF  propagation  only) , 
calculation  flags  for  debris  energy  sources  (determine  whether  ioni- 
zation calculations  will  be  made)  are  determined  and  the  bursts  that 
produce  ionization  above  100  km  are  identified. 

Next,  an  output  option  flag  is  set  to  control  the  amount  of 
detailed  output  and  the  ground  conductivity  for  the  path  set  equal 
to  input  values  are  found  by  calling  routine  CONMAP.  If  neutron  de- 
cay beta  ionization  calculations  are  to  be  made  (input  option),  rou- 
tine EDEPNB  is  initialized.  Then  an  altitude  loop  is  started  and 
the  debris  data  file  and  a data  file  used  for  beta  particle  and  corn- 
ton  electron  data  are  rewound.  The  electron  and  ion  densities  and 
collision  frequencies  at  each  altitude  are  found  by  calling  routine 
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Figure  9.  Flowchart  for  subroutine  ENVIRM. 
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SUBROUTINE  EPRIME 

This  routine  calculates  the  derivatives  with  respect  to  alti- 
tude of  the  anisotropic  reflection  coefficients  for  ELF  propagation. 

The  reflection  coefficients  are  defined  from  a prior  integra- 
tion step.  Use  is  made  of  the  property  that  the  „R|t  and  ARA  values 
for  ELF  cases  approximate  1 and  -1  respectively  (see  Section  2, 

Volume  3).  If  the  earth's  magnetic  field  effects  are  to  be  ignored, 
simple  relationships  defining  the  derivatives  are  used.  More  complex 
definitions  are  used  for  the  case  where  magnetic  field  effects  are 
included. 


SUBROUTINE  ERINT 


This  routine  calculates  ELF  reflection  coefficients  using 
the  Runge-Kutta  integration  algorithms. 

To  avoid  the  possibility  of  interpolation-induced  numerical 
problems,  the  integration  intervals  always  start  and  terminate  at 
adjacent  field-point  altitudes. 

The  initial  reflection  coefficients  are  defined  followed  by 
definition  of  the  lower  altitude  of  the  first  integration  interval. 
The  integration  step  size  is  defined  as  0.5  km,  and  the  derivatives 
of  the  reflection  coefficients  are  calculated. 

Second-  and  fourth-order  Runge-Kutta  integrations  are  per- 
formed. A measure  of  the  error  is  made  by  computing  the  root  mean 
square  of  the  difference  between  the  two  integrations.  This  error 
is  used  to  determine  the  step  size  for  the  next  integration  step, 
which  is  halved,  held  constant,  or  doubled  depending  on  the  magni- 
tude of  the  error. 

When  the  step  size  is  to  remain  constant  or  be  halved,  a 
test  is  made  to  determine  if  the  bottom  of  the  integration  has  been 
reached.  If  it  has,  control  is  returned  to  the  calling  routine.  If 
more  field  altitudes  remain  to  be  processed,  routine  ERINT  does  the 
next  integration  step. 
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SUBROUTINE  ESMAT 

This  routine  calculates  the  S-matrix  of  the  coefficients  for 
the  ELF  reflection  coefficient  differential  equation. 

The  equations  used  for  determining  the  S-matrix  values  are 
in  Appendix  A,  Volume  3. 


SUBROUTINE  ESTART 


This  routine  calculates  the  starting  values  for  ELF  reflection 
coefficient  integration  through  an  anisotropic  ionosphere. 

The  equations  used  to  determine  the  initial  reflection  coeffi- 
cient values  are  defined  in  Appendix  B,  Volume  3.  For  isotropic  ion- 
ospheres the  coupled  terms  tR„  and  MRA  are  zero,  and  nRn  and  ARA  can 
be  determined  from  simple  expressions. 

For  the  anisotropic  case  the  Booker  quartic  coefficients  are 
calculated,  followed  by  solving  for  the  roots  of  the  quartic.  These 
solutions  are  then  used  as  initial  values  for  a Newton-Raphson  itera- 
tion technique  that  refines  the  values.  The  two  roots  that  corres- 
pond to  the  upgoing  waves  are  selected  and  used  to  determine  the  re- 
flection coefficient  matrix  values. 


SUBROUTINE  ETGAIN 

This  routine  calculates  ELF  height  gain  functions. 

Flat-earth  expressions  defined  in  Section  3,  Volume  3,  are 
used  to  calculate  the  ELF  height  gain  functions. 


SUBROUTINE  EXCFAC 


This  routine  calculates  VLF  anisotropic  excitation  factors 

The  excitation  factors  are  determined  as  a function  of  the 
ground  reflection  coefficients  and  eigenangle  from  relations  given 
in  Section  3,  Volume  3. 


SUBROUTINE  FDELH 


This  routine  computes  ionospheric  incidence  angles,  semi-hop 
ground  distance,  and  the  ray  path  length  for  LF  propagation. 

Standard  curved-earth  geometry  formulations  are  used. 


SUBROUTINE  FINDHP 

This  routine  finds  the  reflection  altitude  for  a ray  with  a 
specified  ground  elevation  angle.  A flowchart  for  the  routine  is 
shown  in  Figure  10. 

This  routine  is  called  in  an  iterative  ray-geometry  calcula- 
tion, with  a trial  ground  elevation  angle  defined.  A nominal  reflec- 
tion altitude  and  the  ionospheric  scale  height  and  ionospheric  bottom 
altitude  (below  which  ionization  is  negligible)  are  defined  at  all 
vertical  paths  along  the  great  circle  path. 

Vertical  path  indices  that  bracket  the  reflection  region  are 
defined,  and  a conservative  estimate  of  the  lowest  possible  reflec- 
tion altitude  is  made.  Then  an  iteration  loop  is  initiated  to  find 
the  actual  reflection  height.  Within  the  iteration  loop  the  follow- 
ing steps  occur: 

1.  The  distance  from  the  ray  origin  to  a point  on  the  ray 
that  is  at  the  trial  altitude  is  computed,  given  the 
ray  ground  elevation  angle.  The  ionospheric  incidence 
angle  is  also  computed. 

2.  Interpolation  is  used  to  define  a nominal  reflection 
altitude  and  ionospheric  scale  height. 

3.  The  actual  reflection  altitude  is  computed  using  the 
interpolated  values  from  Step  2 and  the  computed  iono- 
spheric incidence  angle. 

4.  The  trial  altitude  is  incremented  upward  and  Steps  1 
through  4 repeated  until  the  trial  altitude  differs  by 
a specified  minimum  value  from  the  computed  reflection 
altitude. 

5.  The  value  computed  for  the  reflection  altitude  in  the 
last  iteration  is  defined  to  be  the  required  reflection 
altitude. 
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Figure  10.  Flowchart  for  subroutine  FINDHP. 


The  nominal  reflection  altitude  and  ionospheric  scale  height 
and  the  definition  of  the  reflection  altitude  in  terms  of  these  ijuan- 
tities  and  the  ionospheric  incidence  angle  are  explained  and  formulas 
given  in  Section  4,  Volume  3. 


SUBROUTINE  FINGEO 


This  routine  computes  hop  geometry  between  transmitter  and 
receiver  antennas.  A flowchart  for  the  routine  is  shown  in  Figure  11. 

Subroutine  FINGEO  contains  iteration  loops  for  defining  sky- 
wave  geometry  when  the  ionospheric  reflection  altitude  varies  along 
the  propagation  path.  Five  skywave  paths  are  sought  with  N - 2, 

N - 1,  N,  N + 1,  and  N + 2 ionospheric  reflections,  where  N is  the 
minimum  number  of  hops  required  to  traverse  the  total  distance  be- 
tween transmitter  and  receiver  without  diffraction.  When  N - 2 or 
N - 1 is  zero  or  negative,  the  corresponding  path  is  not  considered. 

A maximum  limit  on  the  number  of  hops  is  set  at  nine. 

The  iteration  starts  with  a one-hop  path.  The  ground  eleva- 
tion angle  is  initialized  to  zero,  and  the  total  distance  traveled 
by  a one-hop  ray  is  computed  in  FINDHP.  If  the  distance  traveled  is 
less  than  the  path  distance,  the  hop  counter  is  incremented  and  the 
path  geometry  for  a two-hop  path  is  sought.  This  process  is  continued 
until  the  sum  of  the  hop  distance  exceeds  the  total  path  distance. 

Then  the  ground  elevation  angle  is  incremented  in  an  iteration  loop 
to  make  the  total  hop  distance  agree  with  the  total  path  length  with- 
in a specified  tolerance. 

The  number  of  hops  that  produces  the  first  hop- length  sum 
that  exceeds  the  path  distance  defines  N (described  above).  The  maxi- 
mum number  of  hops  is  defined  as  N + 2,  and  the  diffracted  path  with 
the  greatest  number  of  hops  is  defined  as  N - 1. 

After  defining  N,  the  hop  geometry  for  N + 1 and  N + 2 hop 
skywaves  is  defined  by  iterating  on  ground  elevation  angle. 

Finally,  after  defining  the  geometry  for  all  the  rays  where 
geometric  optics  apply,  the  geometry  for  diffracted  paths  is  computed 
when  N exceeds  1.  The  diffraction  path  geometry  is  obtained  by  di- 
viding the  path  into  equal  increments,  assuming  the  ionospheric  re- 
flection point  occurs  at  the  middle  of  the  increment  and  that  the 
ground  elevation  angle  is  zero. 


Figure  11.  Flowchart  for  subroutine  FINGEO. 


SUBROUTINE  FMAG 


This  routine  computes  horizontal  and  vertical  groundwave  loss 
functions,  and  diffraction  losses.  A flowchart  for  the  routine  is 
shown  in  Figure  12. 

Subroutine  FMAG  computes  all  parameters  used  to  define  the 
influence  of  the  ground  on  the  skywave  fields.  The  groundwave  loss 
term  is  computed  using  Airy  functions  representations.  The  antenna 
foreground  factors  are  computed  using  either  a geometric  representa- 
tion for  ground  elevation  angles  greater  than  2 degrees  or  Airy  func- 
tion formulations  for  smaller  ground  elevation  angles.  Ground  reflec- 
tion loss  or  diffraction  loss  is  computed  at  intermediate  ground  re- 
flection points. 

Both  geometrical  optical  approximations  (Fresnel  coefficients) 
and  the  Airy  function  formulations  are  described  in  detail  in  Section 
4,  Volume  3. 


Figure  12.  Flowchart  for  subroutine  FMAG. 
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SUBROUTINE  GNDWV 


This  routine  computes  the  groundwave  attenuation  function  for 
a specified  path  length,  wave  frequency,  and  ground  impedance. 

An  average  value  for  the  ground  conductivity  is  first  computed. 
The  groundwave  loss  function  for  a homogeneous  curved  earth  is  obtained 
from  Subroutine  FMAG.  The  formulas  evaluated  in  FMAG  are  described  in 
Section  4 of  Volume  3. 
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SUBROUTINE  GREFL 

!;  i 

This  routine  computes  the  Fresnel  ground  reflection  coeffi- 
cient. 

This  routine  computes  the  ground  reflection  coefficient  using 
standard  Fresnel  formulas  for  both  horizontal  and  vertical  polariza- 
tion. It  is  used  for  ground  elevation  angles  greater  than  2 degrees. 
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SUBROUTINE  GUESS 

This  routine  determines  isotropic  mode  solutions  for  VLF 
propagation. 

First,  the  polarization  (horizontal  or  vertical)  is  deter- 
mined and  routine  ISOREF  is  called  to  determine  the  isotropic  reflec- 
tion coefficients.  Then,  a call  to  routine  APPROX  is  made  to  deter- 
mine the  eigenangle.  The  eigenangle  is  referenced  to  the  altitude 
where  the  index  of  refraction  is  assumed  to  be  unity  and  excitation 
and  height  gain  factors  are  determined. 


SUBROUTINE  HTGAIN 


This  routine  calculates  the  anisotropic  antenna  height-gain 
factors  for  VLF  propagation. 

The  VLF  height  gain  factors  are  calculated  as  defined  by 
equations  in  Section  3,  Volume  3.  The  first  time  routine  HTGAIN  is 
called,  an  initialization  is  performed  for  all  parameters  independent 
of  antenna  altitude.  For  subsequent  calls  to  routine  HTGAIN,  only 
altitude-dependent  parameters  are  computed.  If  the  imaginary  part 
of  the  eigenvalue  that  is  input  to  HTGAIN  exceeds  10  degrees,  flat- 
earth  approximations  are  used.  Otherwise  curved-earth  calculations 
are  made. 


SUBROUTINE  INPUT 


This  routine  reads  input  quantities  for  the  WIiDCOM  IV  code, 
converts  input  problem  geometry  into  coordinate  systems  used  within 
the  code,  performs  preliminary  calculations,  and  writes  out  a descrip- 
tion of  the  problem  to  be  run. 

A detailed  description  of  the  input  data  format  and  input 
quantities  is  given  in  Section  2 of  Volume  1.  First,  input  data  (Input 
Blocks  1 through  9)  are  read.  If  input  data  for  particular  Input 
Blocks  are  not  given,  previously  specified  data  (stacked  problem 
cases)  or  default  data  are  used. 

In  conjunction  with  reading  input  data,  the  number  of  calcula- 
tion times  are  determined  if  these  inputs  are  specified.  Also,  if 
weapon  spectral  data  are  given,  routines  W0G1,  W0N1,  and  W0X1  are 
called  to  process  the  data  for  use  in  the  WEDCOM  code. 

After  the  input  data  are  read  in,  a printout  describing  the 
problem  case  is  prepared  and  geometry  calculations  for  the  transmitter 
and  receiver  terminal  locations  and  burst  point  locations  determined. 
The  geometry  calculations  are  performed  to  convert  the  input  geometry 
specification  to  an  earth-centered  geographic  vector  description  and 
to  provide  output  in  all  of  the  several  allowable  input  geometry 
systems . 

If  in  a stacked  problem  case  the  location  of  the  origin  is 
changed  and  the  location  of  the  transmitter  and  receiver  terminal  are 
not  (input  group  6 changed  and  input  group  7 not  given) , the  terminal 
locations  will  be  taken  as  the  geographic  locations  determined  in 
the  problem  case  they  were  specified.  Thus,  even  if  the  terminal  loca- 


tions were  specified  in  relation  to  the  origin  they  will  not  move 
unless  the  input  cards  (input  group  7)  are  given.  The  same  comment 
applies  to  burst  locations. 
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SUBROUTINE  IONOSU 


This  routine  computes  properties  of  the  ambient  ionosphere. 

The  routine  is  a modified  version  of  routine  IONOSU  developed  by 
Science  Applications  Incorporated  for  ROSCOE  (Reference  6) . 

Inputs  to  the  routine  are  the  altitude  of  interest,  time  of 
day  (day  or  night)  , and  atmospheric  properties  (neutral  species  con- 
centrations and  gas  temperature)  at  the  altitude  of  interest.  For 

altitudes  below  90  km  the  only  output  quantity  computed  is  the  ambient 

-3  -1 

ion-pair  production  rate  (cm  sec  ).  For  higher  altitudes,  the 
ambient  electron  density,  atomic  ion  density,  molecular  ion  density, 
and  electron  temperature  are  computed  in  addition  to  the  ambient  ion- 
pair  production  rate.  The  output  quantities  arc  determined  by  fits  to 
nominal  mid-latitude  daytime  and  nighttime  data  and  are  not  a function 
of  geographic  location,  or  date. 

For  use  in  WEDCOM  the  nighttime  electron  density  profile  has 
been  modified  to  have  a minimum  value  between  the  F region  and  the  E 
region  maximum  similar  to  that  in  the  nighttime  profile  adapted  by 
the  Naval  Ocean  Systems  Center  (Reference  4) . 
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SUBROUTINE  ISOREF 

This  routine  computes  coefficients  for  an  exponential  repre- 
sentation of  the  VLF  reflection  coefficient  as  a function  of  the  angle 
of  incidence  and  a phase  reference  altitude  for  a specified  vertical 
ionization  profile.  A flowchart  for  the  routine  is  shown  in  Figure  13. 

A reference  altitude  is  defined  as  the  altitude  where  reflec- 
tion maximizes.  Starting  with  the  lowest  field-point  altitude,  compu- 
tations are  made  of  the  real  and  imaginary  parts  of  the  square  of  the 
complex  index  of  refraction  (Section  3,  Volume  3).  An  exponential 
interpolation  at  1-km  intervals  is  made  between  field-point  altitudes. 
The  reference  altitude  is  defined  as  the  altitude  where  the  imaginary 
part  of  the  square  of  the  index  of  refraction  is  equal  to  0.04.  The 
highest  and  lowest  altitudes  on  the  ionization  profile  that  must  be 
considered  are  defined  prior  to  the  call  to  ISOREF. 

After  the  reference  altitude  is  determined,  the  reflection 
coefficient  for  two  specified  angles  of  incidence  are  found  by  a 
recursive  relationship  that  proceeds  downward  from  the  highest  altitude. 
"Slab"  impedance  values  are  computed  for  1-km  slabs  by  interpolating 
the  index  of  refraction  between  field  points.  Depending  on  the  value 
of  KPOL,  (an  input  option  indicating  polarization)  the  slab  impedances 
are  calculated  for  either  horizontal  or  vertical  polarization.  Com- 
putations of  the  phase  of  the  reflection  coefficients  are  made  at 
altitudes  below  the  reference  altitude. 

The  coefficients  used  in  the  exponential  representation  of 
the  reflection  coefficient  as  a function  of  the  angle  of  incidence 
are  then  computed.  Generally,  reflection  coefficient  calculations 
are  terminated  below  an  altitude  where  significant  absorption  occurs. 
However,  there  is  the  possibility  that  the  ionization  below  the 
reference  altitude  may  decrease  sufficiently  to  cause  the  slab  calcu- 
lation to  be  terminated  and  then  increase  enough  (e.g.,  in  a low-alti- 
tude debris  region  where  the  ionization  does  not  cause  reflection)  to 
produce  some  absorption.  To  allow  for  this  possibility,  absorption 
occuring  below  the  termination  altitude  is  used  to  modify  the  reflection 
coefficient  amplitude. 
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COMPUTE  COEFFICIENTS  FOR 
ANALYTICAL  REPRESENTATION 
OF  REFLECTION  COEFFICIENT 


RETURN 


Figure  13.  Flowchart  for  subroutine  ISOREF. 
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SUBROUTINE  ISORLF 

This  routine  computes  the  amplitude  and  phase  of  the  LF  iso- 
tropic reflection  coefficient  and  a phase  reference  altitude  for  a 
specified  vertical  ionization  profile  and  angle  of  incidence. 

Routine  ISORLF  is  essentially  the  same  as  routine  ISOREF 
except  that  calculations  are  made  for  a specified  angle  of  incidence 
given  as  input  and  the  reflection  coefficient  is  given  as  output 
rather  than  the  coefficients  used  in  an  analytical  representation  of 
the  reflection  coefficient. 
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SUBROUTINE  LENTRP 


This  routine  interpolates  linearly  between  data  stored  for 
adjacent  vertical  paths  to  obtain  electron  and  positive  ion  density 
values,  ground  electrical  constants,  and  geometric  parameters. 

Values  for  ground  conductivity,  ground  relative  dielectric 
constant,  propagation  azimuth,  magnetic  field  strength,  magnetic  field 
dip  angle,  ionospheric  top  and  bottom  altitudes,  ionosphere  scale 
height,  and  electron  and  positive  ion  density  vertical  profiles  are 
stored  for  each  vertical  path.  For  a given  position  between  two 
vertical  paths,  output  quantities  are  found  by  linearly  interpolating 
with  distance. 


SUBROUTINE  LFHOP 


i 

i! 
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This  routine  is  the  LF  propagation  driver  routine  and  deter- 
mines anisotropic  reflection  coefficients  and  ground  wave  and  sky- 
wave  field  strengths.  A simplified  flowchart  for  the  routine  is 
shown  in  Figure  14. 

First,  propagation  parameters  and  geometry  are  defined  and 
nominal  system  quantities  are  written  out.  Then,  for  propagation  paths 
less  than  6000  km  the  groundwave  field  strength  is  determined  by 
calling  routine  GNDWV. 

A time  loop  over  the  calculation  times  is  started  and  an 
ambient  calculation  flag  (determine  whether  calculations  will  be 
made  for  ambient  conditions)  is  set  from  input  parameters.  Next,  a 
loop  over  the  vertical  paths  between  transmitter  and  receiver  is 
started.  For  each  path  magnetic  field  quantities  are  established  and 
the  level  of  detailed  output  requested  determined.  The  ionization 
and  collision  frequencies  along  the  vertical  path  for  disturbed  con- 
ditions are  obtained  from  the  environment  file.  If  detailed  output 
concerning  the  index  of  refraction  has  been  requested  (10UTPG  = 2), 
routine  OUTION  is  called. 

After  the  ionospheric  conditions  for  each  vertical  path 
have  been  obtained  the  skywave  fields  are  determined  for  ambient 
(if  requested)  and  disturbed  conditions.  First,  routine  BHRHS  is 
called  to  obtain  starting  and  stopping  altitudes  for  reflection  coef- 
ficient calculations  and  routine  FINEGO  is  called  to  determine  sky- 
wave  geometry.  Then,  routine  EHOP  is  called  to  determine  the  sky- 
wave  field  strengths.  The  total  field  (sum  of  groundwave  field  and 
skywave  fields)  is  determined  by  calling  routine  TFIELD.  The  last 
step  before  starting  calculations  for  the  next  calculation  time  is 
to  prepare  nominal  output  for  disturbed  and  ambient  field  strengths. 
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Figure  14.  Flowchart  for  subroutine  LFHOP. 
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SUBROUTINE  LINKER 


This  routine  is  used  when  operating  the  WEDCOM  IV  code  on  a 
Honeywell  computer  system  to  provide  overlaying.  The  LINKER  routine 
is  used  rather  than  direct  calls  to  LLINK  (a  Honeywell  linking 
routine)  in  order  to  be  able  to  test  whether  the  overlay  is  already 
available. 


SUBROUTINE  MATINV 


This  routine  solves  a set  of  linear  simultaneous  equations 
[A]X  = B 

where 

[A]  = a square  matrix 
B = force  element  vector 

X = independent  variable  element  vector  to  be  solved. 


Craut's  L - U decomposition  algorithm  is  used  with  partial  pivoting 
(Reference  7) . 


SUBROUTINE  MDHNKL 


This  routine  computes  the  modified  Hankel  functions  of  order 
one-third  and  their  derivatives  for  a specified  complex  argument. 
This  is  an  NOSC  routine  adapted  for  WtDCOM  (Reference  5). 
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SUBROUTINE  MMATRX 

This  routine  calculates  the  ionospheric  susceptibility  matrix 
including  magnetic  field  effects  and  earth  curvature. 


The  squared  ratios  of  the  electron,  positive  ion,  and  negative 
ion  plasma  frequencies  to  the  signal  frequency  are  defined.  Next, 
the  corresponding  ratios  of  collision  frequencies  to  the  signal  fre- 
quency are  defined.  The  susceptibility  matrix  is  initialized. 

The  susceptibility  matrix  is  calculated  by  separately  calcu- 
lating and  then  summing  the  contributions  due  to  electrons  and  posi- 
tive and  negative  ions.  If  field  effects  are  not  important,  only 
the  diagonal  terms  are  calculated.  The  effects  of  earth  curvature 
are  included  by  an  artificial  modification  of  the  diagonal  terms  in 
the  susceptibility  matrix. 
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SUBROUTINE  MODCON 

This  routine  computes  the  VLF  model  mode  conversion  coeffi- 
cients used  to  determine  the  received  signal  strengths  from  non-uniform 
waveguides . 

A test  is  made  to  determine  if  the  waveguide  slab  characteris- 
tics for  the  slab  under  consideration  are  identical  to  the  prior  slab. 

If  so,  a simple  relationship  exists  between  the  coefficients  of  the 
new  slab  and  the  prior  slab.  For  the  first  slab,  or  for  the  case  where 
two  adjacent  slab  characteristics  differ,  height-gain  characteristics 
(in  Volume  3 Equations  3-65  through  3-70)  at  the  slab  interface  are 
computed  using  the  eigenangles  for  the  new  slab.  Slab  interface  height- 
gain  characteristics  are  also  computed  for  the  prior  slab  for  the  same 
eigenangles.  These  characteristics  are  then  used  to  formulate  boundary 
condition  relationships  (Equation  3-64) . The  mode  conversion  coeffi- 
cients are  computed  from  the  solution  of  a set  of  simultaneous  equations 
(Equation  3-63) . The  height-gain  characteristics  computed  above  and 
mode  conversion  coefficients  are  saved  for  the  next  slab  interface. 
Detailed  output  is  printed  if  desired. 
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SUBROUTINE  MODEZ 

This  routine  computes  solutions  to  the  VLF  mode  equation  for 
a given  ionospheric  profile  and  ground  conductivity.  Height  gain  and 
excitation  factors  are  also  calculated.  A simplified  flowchart  for 
the  routine  is  shown  in  Figure  14. 

For  ground  ranges  less  than  4000  km  and  frequencies  less  than 
15  khz  six  moucs  are  computed;  otherwise  four  are  computed.  A flag 
is  initialized  at  the  start  of  the  mode  loop  to  direct  calculations 
to  the  TM  or  TE  mode  branch,  half  of  the  mode  solutions  result  from 
I'M  calculations  and  the  remaining  half  from  TE  calculations.  An 
isotropic  solution  is  found  by  calling  routine  GUESS  and  the  attenua- 
tion rate  is  computed.  If  the  attenuation  rate  is  greater  than  20  dB/Mm, 
a logical  parameter,  ISOANS,  is  set  to  1 to  indicate  that  the  isotropic 
solution  will  be  used. 

Routine  EIGENV  is  employed  to  calculate  the  value  of  the  eigen- 
value, including  earth  curvature  and  the  earth's  magnetic  field  effects, 
using  the  approximate  isotropic  ionosphere  eigenvalue  as  a first  estimate 
in  the  iteration  process.  If  the  iteration  process  does  not  converge, 
the  isotropic  value  is  saved  for  further  calculations  and  a new  mode 
is  analyzed. 

If  the  logical  parameter  ISOANS  is  equal  to  1 (either  because 
of  the  attenuation  rate  or  because  of  being  set  in  routine  VLFMOD) , 
the  isotropic  solution  obtained  from  routine  GUESS  is  used  and  only 
parameters  related  to  mode  conversion  are  obtained  from  routine  EIGENV. 

After  the  eigenvalue  is  calculated,  the  height  gain  and  exci- 
tation factors  are  found  by  calling  routines  HTGAIN  and  EXCFAC.  After 
all  the  modes  are  calculated  the  mode  outputs  are  reordered  so  that 
the  real  part  of  the  eigenangles  decreases  as  the  mode  number  increases. 

If  detailed  output  has  been  requested,  the  mode  quantities, 
excitation  factors  and  height-gain  factors  are  written  out. 
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Figure  15.  Flowchart  for  subroutine  MODEZ. 


SUBROUTINE  OUTION 


This  routine  computes  the  ionospheric  index  of  refraction 
profile,  the  reference  altitude  for  VLF  and  LF  reflection,  and  the  top 
and  bottom  altitudes  for  reflection  coefficient  calculations.  A flow- 
chart for  the  routine  is  shown  in  Figure  16. 

The  anisotropic  and  isotropic  indexes  of  refraction  are  defined 
at  each  field  point  beginning  at  the  lowest  field  point.  The  aniso- 
tropic index  is  used  for  detailed  output  and  determining  the  top  alti- 
tude for  ELF  reflection  coefficient  calculations.  The  isotropic  index 
is  used  for  approximating  the  VLF  and  LF  reference  altitudes  and  the 
top  and  bottom  reflection  coefficient  calculation  altitudes. 


For  ELF  the  top  starting  altitude  is  determined  as  the  altitude 
where  the  one-way  vertical  absorption  equals  30  dB  and  the  bottom  alti- 
tude is  set  to  zero.  For  VL  and  LF  the  top  starting  altitude  is  deter- 
mined as  the  altitude  where 


B = imaginary  part  of  the  square  of  the  index  of  refraction. 
The  bottom  altitude  for  VLF  and  reflection  coefficient  calcu- 
lations is  estimated  as  the  altitude  below  where 
B < 0.0002 

and  the  reference  altitude  is  estimated  as  the  altitudes  where 
B = 0.04  . 


After  completing  the  field-point  altitude  loop,  refinement  of 
the  ionospheric  top  and  reflection  altitude  definitions  is  made  by 
interpolating  between  field  points. 
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Figure  16.  Flowchart  for  subroutine  OUTION. 
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SUBROUTINE  QSUBI 

This  routine  computes  the  normalized  ionospheric  surface  impe- 
dance. 

This  routine  employs  equations  given  in  Section  3,  Volume  4 
to  determine  the  normalized  ionospheric  surface  impedance  assuming 
vertical  or  horizontal  polarization. 


SUBROUTINE  REFLCO 


I 


This  routine  computes  complex  values  of  reflection  coefficients. 
A flowchart  for  the  routine  is  shown  in  Figure  17. 

Subroutine  REFLCO  is  the  driver  routine  for  calculating  iono- 
spheric reflection  coefficients.  Magnetic  field  quantities,  ionospheric 
incident  angle  values,  and  a starting  altitude  for  the  reflection  co- 
efficient integration  are  defined.  Electron  and  ion  densities,  col- 
lision frequencies,  and  elements  of  the  susceptibility  matrix  are 
defined  at  the  starting  altitude  by  calling  routines  ENNU  and  MMATRX. 
Starting  values  for  the  reflection  coefficients  are  obtained  from 
routine  START,  and  the  integration  process  is  initiated  by  calling 
routine  RINT. 


\ 
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SUBROUTINE  REFLCO 


Figure  17.  Flowchart  for  subroutine  REFLCO. 
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SUBROUTINE  REORDER 

This  routine  reorders  environmental  data  stored  in  a file  for 
each  vertical  path  for  all  calculation  times  so  that  they  are  stored 
for  each  calculation  time  for  all  vertical  paths. 


SUBROUTINE  RGND 


— 


This  routine  calculates  ELF  and  VLF  ground  reflection  coeffi- 
cients. 

After  initialization  of  several  conductivity  parameters,  a 
test  determines  if  this  is  a VLF  or  ELF  case.  For  VLF  the  reflection 
coefficients  arc  calculated  using  equations  given  in  Appendix  C, 

Volume  3.  For  ELF  the  Fresnel  ground  reflection  coefficients  are  cal- 
culated. 
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SUBROUTINE  RINT 


This  routine  calculates  reflection  coefficients  using  the 
Runge-Kutta  integration  algorithms.  A flowchart  for  the  routine  is 
shown  in  Figure  18. 

To  avoid  the  possibility  of  interpolation-induced  numerical 
problems,  the  integration  intervals  always  start  and  terminate  at 
adjacent  field-point  altitudes. 

The  reflection  coefficient  starting  values  are  defined  followed 
by  definition  of  the  lower  altitude  of  the  first  integration  interval. 
The  integration  step  size  is  defined  as  0.5  km  and  the  derivatives  of 
the  reflection  coefficients  are  calculated. 

Second-  and  fourth-order  Runge-Kutta  integrations  are  performed. 
A measure  of  the  error  is  made  through  a root-mean  square  of  the  dif- 
ference between  the  two  integrations.  This  error  is  used  to  determine 
the  step  size  for  the  next  integration  step,  which  is  halved,  held 
constant,  or  doubled  depending  on  the  magnitude  of  the  error. 


When  the  step  size  is  to  remain  constant  or  be  halved,  a test 
is  made  to  determine  if  the  bottom  of  the  integration  interval  has 
been  reached.  If  it  has,  control  is  returned  to  the  calling  routine. 
Otherwise,  if  more  field  altitudes  remain  to  be  processed.  Subroutine 
RINT  does  the  next  integration  step. 


Figure  18.  Flowchart  for  subroutine  RINT. 
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SUBROUTINE  ROOT 

This  routine  employs  Newton's  procedure  to  compute  the  root 
of  the  VLF  mode  equation 

Subroutine  ROOT  employs  Newton's  iteration  method  to  calculate 
the  solution  to  the  isotropic  mode  equation  from  a given  approximate 
solution.  The  parameter  AMIN  is  defined  as  the  absolute  part  of 
(1  - AB) , where  A and  B are  defined  in  Section  3,  Volume  7.  To  obtain 
a good  solution,  AMIN  must  be  smaller  than  0.CU5,  and  the  magnitude 
of  the  difference  between  successive  calculations  of  the  root  must  not 
be  larger  than  0.005.  When  these  conditions  are  not  satisfied  after 
20  iterations,  a statement  is  written  indicating  the  last  value  of 
AMIN.  When  the  process  diverges,  the  root  is  set  to  the  original 
guess  if  it  is  a good  approximation;  otherwise  a large  attenuation  is 
assigned. 
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SUBROUTINE  RPRIME 


i 

i 

This  routine  calculates  the  derivatives  with  respect  to  alti- 
tude of  the  anisotropic  reflection  coefficients. 

The  reflection  coefficients  are  defined  from  a prior  integra- 
tion step,  if  the  earth's  magnetic  field  effects  are  to  be  ignored, 
simple  relationships  defining  the  derivatives  are  used.  More  complex 
definitions  are  used  for  the  case  where  magnetic  field  effects  are 
included  (see  Section  2,  Volume  3). 


SUBROUTINE  RSTART 


This  routine  calculates  the  starting  values  for  reflection 
coefficient  integration  through  anisotropic  ionospheres. 

The  equations  used  to  determine  the  initial  reflection  coeffi- 
cient values  are  defined  in  Appendix  B,  Volume  3.  For  isotropic  iono- 
spheres the  coupled  terms  jR|(  and  nRj_  are  zero,  and  i(Rh  and  jR^ 
can  be  determined  from  simple  expressions. 

For  the  anisotropic  case,  the  Booker  quartic  coefficients  are 
calculated,  followed  by  solving  for  the  roots  of  the  quartic.  These 
solutions  arc  then  used  as  initial  values  for  a Newton-Raphson  itera- 
tion technique  that  refines  the  values.  The  two  roots  that  corres- 
pond to  the  upgoing  waves  are  selected  and  used  to  determine  initial 
reflection  coefficient  matrix  values. 
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SUBROUTINE  RUNG 

This  routine  perforins  a Runge-Kutta  integration  in  the  ground- 
wave  calculations. 

A standard  fourth-order  Runge-Kutta  integration  procedure  is 

used.  One  of  the  following  two  equations  is  integrated  over  the 

inverval  from  0 to  Q : 
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SUBROUTINE  SMAT 
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This  routine  calculates  the  S-matrix  of  the  coefficients  for 
the  reflection  coefficient  differential  equation. 

The  equations  used  for  determining  the  S-matrix  values  are  in 
Appendix  A,  Volume  3. 
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SUBROUTINE  SURFQ 

1 his  routine  determines  the  normalized  surface  impedance  for 
horizontally  or  vertically  polarized  fields. 

Subroutine  SURl-'Q  solves  equations  given  in  Section  3,  Volume 
3 to  determine  the  normalized  ground  impedance  for  either  the  horizon- 
tally or  vertically  polarized  fields,  depending  on  the  option  selected 
at  input. 
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SUBROUTINE  TEM 


1 
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This  routine  calculates  TEM  mode  field  values  at  the  receiver 
for  ELF  propagation.  A simplified  flowchart  for  the  routine  is  shown 
in  Figure  19. 

First,  propagation  parameters  and  geometry  are  defined  and  nominal 
system  output  quantities  are  written  out.  Then  a time  loop  over  the 
calculation  times  is  started  and  eigenvalue-dependent  data  parameters 
defined  in  routine  ELFMOD  are  read  from  a data  file.  If  ambient  calcu- 
lations have  been  requested  (input  option)  data  for  ambient  conditions  are 
read  in  first;  otherwise,  data  for  disturbed  conditions  are  obtained. 

The  short-  and  long-path  height-gain  factors  for  the  transmitter  and 
receiver  locations  are  determined  by  calls  to  routine  ETGAIN.  The 
short  path  WKB  approximation  is  always  calculated  and  the  long  path 
calculation  is  made  when  requested  (input  option).  Then,  the  short 
path  (and  the  long  path  when  requested)  electromagnetic  fields  at  the 
receiver  location  are  determined  and  nominal  output  is  written  out. 

If  the  multiple  receiver  location  option  (input  option)  has  been 


exercised,  the  electromagnetic  fields  at  each  vertical  path  location 
are  computed.  After  calculations  have  been  completed  for  disturbed 
conditions,  calculations  are  started  for  the  next  calculation  time. 


CALCULATIONS^ 
COMPLETED  FOR 
DISTURBED 
. CONDITIONS?. 


NEXT  TIME  1 YES 


Figure  19.  Flowchart  for  subroutine  TEM. 
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SUBROUTINE  TFIELD 

This  routine  computes  the  total  LF  electric  field  strength 
by  summing  the  fields  of  all  the  skywaves  and  the  groundwave. 

The  electric  field  is  first  set  equal  to  the  groundwave  value. 
Then,  a vector  or  rms  sum  of  the  groundwave  and  skywave  fields  is 
computed.  The  vector  sum  is  used  when  the  multiple  receiver  location 
option  (input  option)  is  exercised  and  calculations  of  the  electric 
field  are  made  at  each  vertical  path  location.  The  rms  sum  is  used 
when  the  multiple  receiver  location  option  is  not  exercised  and  calcu- 
lations of  the  electric  field  are  only  made  at  one  location.  If 
detailed  output  has  been  requested  (input  option),  the  field  strengths 
for  vertical,  normal,  and  parallel  field  polarizations  at  the  trans- 
mitter and  receiver  locations  are  written  out. 
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SUBROUTINE  UANTD 


This  routine  determines  the  orientation  of  the  transmitter 
antenna  as  a function  of  time. 

This  is  a dummy  routine  for  illustration.  The  transmitter 
antenna  zenith  angle  and  azimuth  angle  are  both  set  to  zero  (vertical 
dipole) . 
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SUBROUTINE  VLFMCV 

This  routine  computes  the  VLF  electric  field  strength  in  each 
mode  at  each  vertical  path  location  and  the  total  field  strength  at 
receiver  locations  using  mode  conversion  to  account  for  effects  of  a 
variable  height  earth-ionosphere  waveguide. 

Parameters  are  first  defined  which  relate  the  transmitter 
power  to  antenna  orientation.  Then  cumulative  mode  conversion  coeffi- 
cients are  computed  for  the  region  between  the  transmitter  (first 
vertical  path  location)  and  the  current  vertical  path  location.  The 
eigenangles,  height-gain,  and  excitation  factors  for  the  first  verti- 
cal path  location  are  saved.  Then,  if  the  current  vertical  path  loca- 
tion is  a receiver  location,  the  total  field  is  computed  from  either 
a vector  or  rms  sum  of  the  mode  fields.  The  vector  sum  is  used  when 
the  multiple  receiver  location  option  is  exercised  (input  option)  and 
calculations  of  the  electric  field  are  made  at  each  vertical  path 
location.  The  rms  sum  is  used  when  only  one  receiver  location  is 
specified. 
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SUBROUTINE  VLFMOD 


I 

This  routine  is  the  VLF  propagation  driver  routine  and  deter- 
mines anisotropic  reflection  coefficients,  mode  solutions,  and  excita- 
tion factors.  A simplified  flowchart  for  the  routine  is  shown  in 
Figure  20. 

First,  propagation  parameters  and  geometry  are  defined  and 
nominal  system  output  quantities  are  written  out.  Then,  a time  loop 
over  the  calculation  times  is  started  and  an  ambient  calculation  flag 
(.determines  whether  calculations  will  be  made  for  ambient  conditions) 
is  set  from  input  parameters. 

Next,  a loop  over  the  vertical  paths  between  transmitter  and 
receiver  is  started  and  calculations  made  for  the  disturbed  environ- 
ment. For  each  path  the  ionization  and  collision  frequencies  are 
obtained  from  the  environment  file.  Magnetic  field  and  ground  conduc- 
tivity quantities  are  established  and  the  level  of  detailed  output 
requested  is  determined.  Then,  routine  0UTI0N  is  called  to  determine 
a reference  altitude  and  the  top  and  bottom  altitudes  to  be  used  in 
computing  reflection  coefficients.  If  the  reference  altitude  is  less 
than  55  km,  a logical  parameter,  1S0ANS,  is  set  to  1 to  indicate  that 
an  isotropic  mode  solution  is  to  be  used. 

If  the  input  option  to  minimize  the  number  of  reflection  coef- 
ficient calculations  is  exercised  (MRCOFF  = 1),  tests  are  made  to 
sec  how  much  the  top  altitude  obtained  from  routine  01JTI0N,  the  magne- 
tic field  properties,  and  the  ground  conductivity  have  changed  from 
the  values  computed  for  the  previous  vertical  path.  If  they  have  not 
changed  by  more  than  specified  amounts,  the  mode  characteristics 
previously  determined  are  used  for  the  current  vertical  path.  Other- 
wise, either  routine  EIC.SCH  (reference  altitudes  greater  than  70  km) 
or  routine  MODFZ  (reference  altitude  less  than  70  km)  is  called  to 
determine  the  mode  characteristics. 


I 
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The  field  strength  of  each  mode  at  the  vertical  path  location 
is  determined  by  either  calling  routine  VLFWKB,  if  the  WKB  approxima- 
tions are  used,  or  calling  routine  VLFMCV,  if  mode  conversion  is  to 
be  determined.  (Either  routine  VLl-'WKB  or  routine  VLFMCV  must  be 
chosen  as  overlay  [7,3]  before  running  WEDCOM  IV).  Calculations  are 
made  for  each  vertical  path  location  for  disturbed  conditions  and 
then  repeated  for  ambient  conditions  if  ambient  calculations  have 
been  requested  for  the  specified  calculation  time. 
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SUBROUTINE  VLFWKB 


This  routine  computes  the  VLF  electric  field  strength  in  each 
mode  (WKB  approximationj  at  each  vertical  path  location  and  the  total 
field  strength  at  receiver  locations. 


First,  the  electric  field  strength  in  each  mode  at  the  current 
vertical  path  location  is  computed  using  a WKB  approximation.  Then, 
if  the  vertical  path  location  is  a receiver  location,  the  total  field 
is  computed  from  either  a vector  or  rms  sum  of  the  mode  fields.  The 
vector  sum  is  used  when  the  multiple  receiver  location  option  is 
exercised  (input  option)  and  calculations  of  the  electric  field  are 
made  at  each  vertical  path  location.  The  rms  sum  is  used  when  only 
one  receiver  location  is  specified. 
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SUBROUTINE  WFCTVL 

This  is  a modified  Naval  Ocean  Systems  Center  routine  used  in 
mode  search  calculations  (Reference  5).  The  routine  computes  F 
function  values  at  mesh  points. 


fi 


SUBROUTINE  WFINDF 


, 


This  is  a modified  NOSC  routine  used  in  the  mode  search  algo- 
rithm. It  is  used  to  compute  the  modified  mode  equation  (Equation 
3-8  in  Volume  3)  at  a corner  of  a particular  mesh. 
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SUBROUTINE  WFDFDT 


Ihis  is  a modified  Naval  Ocean  Systems  Center  routine  used 
in  mode  search  calculations  (Reference  5j.  The  routine  computes  F 
function  values  and  their  derivatives  with  respect  to  6 at  arbitrary 
values  of  0. 
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SUBROUTINE  WFSINT 

This  is  a Naval  Ocean  Systems  Center  routine  used  in  mode 
search  calculations  (Reference  5) . The  routine  performs  an  integra- 
tion of  the  differential  equations  for  the  ionosphere  reflection 
matrix  through  a free  space  region  over  a curved  earth. 

The  routine  has  one  additional  entry  point  besides  the  main 


entry  point.  Entry  point  WINIFS  is  used  for  an  integration  that  is 
independent  of  0. 


SUBROUTINE  WFZERO 


This  is  a modified  Naval  Ocean  Systems  Center  routine  used  in 
mode  search  computations  (Reference  5).  The  routine  finds  the  zeros 
of  a complex  function,  F,  which  are  within  a specified  rectangular 
region  of  the  complex  (0)  plane,  provided  the  function  has  no  poles 
in  the  vicinity  of  the  rectangle. 


SUBROUTINE  WLGRNG 


This  is  a modified  Naval  Ocean  Systems  Center  routine  used  in 
mode  search  calculations  (Reference  5).  The  routine  performs  Lagrange 
interpolation  of  reflection  coefficients  in  the  cos  0 plane. 

The  routine  has  three  additional  entry  points  besides  the  main 
entry  point.  Entry  point  WINILG  is  used  to  initialize  the  Lagrange 
interpolation.  Entry  WLGDER  is  used  to  compute  the  derivatives  of  the 
interpolated  values  of  cos  0.  Entry  WSETLG  is  used  for  further  ini- 
tialization of  the  Lagrange  interpolation. 


SUBROUTINE  WNOMES 


Ihis  is  amodified  Naval  Ocean  Systems  Center  routine  used  in 
mode  search  calculations  (Reference  5) . The  routine  finds  exact  (in 
the  sense  of  no  mesh  approximations)  locations  of  zeros  of  the  funct 
F for  which  a complete,  but  approximate,  set  was  found  in  routine 
WFZERO. 
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SUBROUTINE  WQUAD 

This  is  a modified  Naval  Ocean  Systems  Center  routine  used  in 
mode  search  calculations  (Reference  5).  The  routine  finds  the  solu- 
tion for  the  real  roots  of  a quadratic  equation  of  the  form 

2 

ax“  + 2bx  + c = 0 


SUBROUTINE  WRBARS 

This  is  a modified  Naval  Ocean  Systems  Center  routine  used  in 
mode  search  calculations  (Reference  5) . The  routine  computes  values 
of  variables  used  to  form  the  elements  of  the  reflection  coefficient 
matrix  of  an  ELM  wave  from  the  earth's  surface. 

The  routine  has  one  additional  entry  point  besides  the  main 
entry  point.  Entry  point  WINFRB  is  used  for  computation  that  is  inde- 
pendent of  0. 
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SUBROUTINE  WSETRH 

This  is  a modified  Naval  Ocean  Systems  Center  routine  used  in 
mode  search  calculations  (Reference  5).  The  routine  selects  given  ; 

points  in  the  Lagrange  interpolation  and  selects  the  reference  height. 
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